// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include <cstdio>

#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "ClpGubMatrix.hpp"
//#include "ClpGubDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpGubMatrix::ClpGubMatrix()
  : ClpPackedMatrix()
  , sumDualInfeasibilities_(0.0)
  , sumPrimalInfeasibilities_(0.0)
  , sumOfRelaxedDualInfeasibilities_(0.0)
  , sumOfRelaxedPrimalInfeasibilities_(0.0)
  , infeasibilityWeight_(0.0)
  , start_(NULL)
  , end_(NULL)
  , lower_(NULL)
  , upper_(NULL)
  , status_(NULL)
  , saveStatus_(NULL)
  , savedKeyVariable_(NULL)
  , backward_(NULL)
  , backToPivotRow_(NULL)
  , changeCost_(NULL)
  , keyVariable_(NULL)
  , next_(NULL)
  , toIndex_(NULL)
  , fromIndex_(NULL)
  , model_(NULL)
  , numberDualInfeasibilities_(0)
  , numberPrimalInfeasibilities_(0)
  , noCheck_(-1)
  , numberSets_(0)
  , saveNumber_(0)
  , possiblePivotKey_(0)
  , gubSlackIn_(-1)
  , firstGub_(0)
  , lastGub_(0)
  , gubType_(0)
{
  setType(16);
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpGubMatrix::ClpGubMatrix(const ClpGubMatrix &rhs)
  : ClpPackedMatrix(rhs)
{
  numberSets_ = rhs.numberSets_;
  saveNumber_ = rhs.saveNumber_;
  possiblePivotKey_ = rhs.possiblePivotKey_;
  gubSlackIn_ = rhs.gubSlackIn_;
  start_ = ClpCopyOfArray(rhs.start_, numberSets_);
  end_ = ClpCopyOfArray(rhs.end_, numberSets_);
  lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
  upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
  status_ = ClpCopyOfArray(rhs.status_, numberSets_);
  saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
  savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
  int numberColumns = getNumCols();
  backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
  changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
  // find longest set
  int *longest = new int[numberSets_];
  CoinZeroN(longest, numberSets_);
  int j;
  for (j = 0; j < numberColumns; j++) {
    int iSet = backward_[j];
    if (iSet >= 0)
      longest[iSet]++;
  }
  int length = 0;
  for (j = 0; j < numberSets_; j++)
    length = std::max(length, longest[j]);
  next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
  toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
  infeasibilityWeight_ = rhs.infeasibilityWeight_;
  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
  noCheck_ = rhs.noCheck_;
  firstGub_ = rhs.firstGub_;
  lastGub_ = rhs.lastGub_;
  gubType_ = rhs.gubType_;
  model_ = rhs.model_;
}

//-------------------------------------------------------------------
// assign matrix (for space reasons)
//-------------------------------------------------------------------
ClpGubMatrix::ClpGubMatrix(CoinPackedMatrix *rhs)
  : ClpPackedMatrix(rhs)
  , sumDualInfeasibilities_(0.0)
  , sumPrimalInfeasibilities_(0.0)
  , sumOfRelaxedDualInfeasibilities_(0.0)
  , sumOfRelaxedPrimalInfeasibilities_(0.0)
  , infeasibilityWeight_(0.0)
  , start_(NULL)
  , end_(NULL)
  , lower_(NULL)
  , upper_(NULL)
  , status_(NULL)
  , saveStatus_(NULL)
  , savedKeyVariable_(NULL)
  , backward_(NULL)
  , backToPivotRow_(NULL)
  , changeCost_(NULL)
  , keyVariable_(NULL)
  , next_(NULL)
  , toIndex_(NULL)
  , fromIndex_(NULL)
  , model_(NULL)
  , numberDualInfeasibilities_(0)
  , numberPrimalInfeasibilities_(0)
  , noCheck_(-1)
  , numberSets_(0)
  , saveNumber_(0)
  , possiblePivotKey_(0)
  , gubSlackIn_(-1)
  , firstGub_(0)
  , lastGub_(0)
  , gubType_(0)
{
  setType(16);
}

/* This takes over ownership (for space reasons) and is the
   real constructor*/
ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix *matrix, int numberSets,
  const int *start, const int *end,
  const double *lower, const double *upper,
  const unsigned char *status)
  : ClpPackedMatrix(matrix->matrix())
  , sumDualInfeasibilities_(0.0)
  , sumPrimalInfeasibilities_(0.0)
  , sumOfRelaxedDualInfeasibilities_(0.0)
  , sumOfRelaxedPrimalInfeasibilities_(0.0)
  , numberDualInfeasibilities_(0)
  , numberPrimalInfeasibilities_(0)
  , saveNumber_(0)
  , possiblePivotKey_(0)
  , gubSlackIn_(-1)
{
  model_ = NULL;
  numberSets_ = numberSets;
  start_ = ClpCopyOfArray(start, numberSets_);
  end_ = ClpCopyOfArray(end, numberSets_);
  lower_ = ClpCopyOfArray(lower, numberSets_);
  upper_ = ClpCopyOfArray(upper, numberSets_);
  // Check valid and ordered
  int last = -1;
  int numberColumns = matrix_->getNumCols();
  int numberRows = matrix_->getNumRows();
  backward_ = new int[numberColumns];
  backToPivotRow_ = new int[numberColumns];
  changeCost_ = new double[numberRows + numberSets_];
  keyVariable_ = new int[numberSets_];
  // signal to need new ordering
  next_ = NULL;
  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
    backward_[iColumn] = -1;

  int iSet;
  for (iSet = 0; iSet < numberSets_; iSet++) {
    // set key variable as slack
    keyVariable_[iSet] = iSet + numberColumns;
    if (start_[iSet] < 0 || start_[iSet] >= numberColumns)
      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
    if (end_[iSet] < 0 || end_[iSet] > numberColumns)
      throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
    if (end_[iSet] <= start_[iSet])
      throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix");
    if (start_[iSet] < last)
      throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix");
    last = end_[iSet];
    int j;
    for (j = start_[iSet]; j < end_[iSet]; j++)
      backward_[j] = iSet;
  }
  // Find type of gub
  firstGub_ = numberColumns + 1;
  lastGub_ = -1;
  int i;
  for (i = 0; i < numberColumns; i++) {
    if (backward_[i] >= 0) {
      firstGub_ = std::min(firstGub_, i);
      lastGub_ = std::max(lastGub_, i);
    }
  }
  gubType_ = 0;
  // adjust lastGub_
  if (lastGub_ > 0)
    lastGub_++;
  for (i = firstGub_; i < lastGub_; i++) {
    if (backward_[i] < 0) {
      gubType_ = 1;
      printf("interior non gub %d\n", i);
      break;
    }
  }
  if (status) {
    status_ = ClpCopyOfArray(status, numberSets_);
  } else {
    status_ = new unsigned char[numberSets_];
    memset(status_, 0, numberSets_);
    int i;
    for (i = 0; i < numberSets_; i++) {
      // make slack key
      setStatus(i, ClpSimplex::basic);
    }
  }
  saveStatus_ = new unsigned char[numberSets_];
  memset(saveStatus_, 0, numberSets_);
  savedKeyVariable_ = new int[numberSets_];
  memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
  noCheck_ = -1;
  infeasibilityWeight_ = 0.0;
}

ClpGubMatrix::ClpGubMatrix(const CoinPackedMatrix &rhs)
  : ClpPackedMatrix(rhs)
  , sumDualInfeasibilities_(0.0)
  , sumPrimalInfeasibilities_(0.0)
  , sumOfRelaxedDualInfeasibilities_(0.0)
  , sumOfRelaxedPrimalInfeasibilities_(0.0)
  , infeasibilityWeight_(0.0)
  , start_(NULL)
  , end_(NULL)
  , lower_(NULL)
  , upper_(NULL)
  , status_(NULL)
  , saveStatus_(NULL)
  , savedKeyVariable_(NULL)
  , backward_(NULL)
  , backToPivotRow_(NULL)
  , changeCost_(NULL)
  , keyVariable_(NULL)
  , next_(NULL)
  , toIndex_(NULL)
  , fromIndex_(NULL)
  , model_(NULL)
  , numberDualInfeasibilities_(0)
  , numberPrimalInfeasibilities_(0)
  , noCheck_(-1)
  , numberSets_(0)
  , saveNumber_(0)
  , possiblePivotKey_(0)
  , gubSlackIn_(-1)
  , firstGub_(0)
  , lastGub_(0)
  , gubType_(0)
{
  setType(16);
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpGubMatrix::~ClpGubMatrix()
{
  delete[] start_;
  delete[] end_;
  delete[] lower_;
  delete[] upper_;
  delete[] status_;
  delete[] saveStatus_;
  delete[] savedKeyVariable_;
  delete[] backward_;
  delete[] backToPivotRow_;
  delete[] changeCost_;
  delete[] keyVariable_;
  delete[] next_;
  delete[] toIndex_;
  delete[] fromIndex_;
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpGubMatrix &
ClpGubMatrix::operator=(const ClpGubMatrix &rhs)
{
  if (this != &rhs) {
    ClpPackedMatrix::operator=(rhs);
    delete[] start_;
    delete[] end_;
    delete[] lower_;
    delete[] upper_;
    delete[] status_;
    delete[] saveStatus_;
    delete[] savedKeyVariable_;
    delete[] backward_;
    delete[] backToPivotRow_;
    delete[] changeCost_;
    delete[] keyVariable_;
    delete[] next_;
    delete[] toIndex_;
    delete[] fromIndex_;
    numberSets_ = rhs.numberSets_;
    saveNumber_ = rhs.saveNumber_;
    possiblePivotKey_ = rhs.possiblePivotKey_;
    gubSlackIn_ = rhs.gubSlackIn_;
    start_ = ClpCopyOfArray(rhs.start_, numberSets_);
    end_ = ClpCopyOfArray(rhs.end_, numberSets_);
    lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
    upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
    status_ = ClpCopyOfArray(rhs.status_, numberSets_);
    saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
    savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
    int numberColumns = getNumCols();
    backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
    changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
    // find longest set
    int *longest = new int[numberSets_];
    CoinZeroN(longest, numberSets_);
    int j;
    for (j = 0; j < numberColumns; j++) {
      int iSet = backward_[j];
      if (iSet >= 0)
        longest[iSet]++;
    }
    int length = 0;
    for (j = 0; j < numberSets_; j++)
      length = std::max(length, longest[j]);
    next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
    toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
    sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    infeasibilityWeight_ = rhs.infeasibilityWeight_;
    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    noCheck_ = rhs.noCheck_;
    firstGub_ = rhs.firstGub_;
    lastGub_ = rhs.lastGub_;
    gubType_ = rhs.gubType_;
    model_ = rhs.model_;
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase *ClpGubMatrix::clone() const
{
  return new ClpGubMatrix(*this);
}
/* Subset clone (without gaps).  Duplicates are allowed
   and order is as given */
ClpMatrixBase *
ClpGubMatrix::subsetClone(int numberRows, const int *whichRows,
  int numberColumns,
  const int *whichColumns) const
{
  return new ClpGubMatrix(*this, numberRows, whichRows,
    numberColumns, whichColumns);
}
/* Returns a new matrix in reverse order without gaps
   Is allowed to return NULL if doesn't want to have row copy */
ClpMatrixBase *
ClpGubMatrix::reverseOrderedCopy() const
{
  return NULL;
}
int ClpGubMatrix::hiddenRows() const
{
  return numberSets_;
}
/* Subset constructor (without gaps).  Duplicates are allowed
   and order is as given */
ClpGubMatrix::ClpGubMatrix(
  const ClpGubMatrix &rhs,
  int numberRows, const int *whichRows,
  int numberColumns, const int *whichColumns)
  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
{
  // Assuming no gub rows deleted
  // We also assume all sets in same order
  // Get array with backward pointers
  int numberColumnsOld = rhs.matrix_->getNumCols();
  int *array = new int[numberColumnsOld];
  int i;
  for (i = 0; i < numberColumnsOld; i++)
    array[i] = -1;
  for (int iSet = 0; iSet < numberSets_; iSet++) {
    for (int j = start_[iSet]; j < end_[iSet]; j++)
      array[j] = iSet;
  }
  numberSets_ = -1;
  int lastSet = -1;
  bool inSet = false;
  for (i = 0; i < numberColumns; i++) {
    int iColumn = whichColumns[i];
    int iSet = array[iColumn];
    if (iSet < 0) {
      inSet = false;
    } else {
      if (!inSet) {
        // start of new set but check okay
        if (iSet <= lastSet)
          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
        lastSet = iSet;
        numberSets_++;
        start_[numberSets_] = i;
        end_[numberSets_] = i + 1;
        lower_[numberSets_] = lower_[iSet];
        upper_[numberSets_] = upper_[iSet];
        inSet = true;
      } else {
        if (iSet < lastSet) {
          throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
        } else if (iSet == lastSet) {
          end_[numberSets_] = i + 1;
        } else {
          // new set
          lastSet = iSet;
          numberSets_++;
          start_[numberSets_] = i;
          end_[numberSets_] = i + 1;
          lower_[numberSets_] = lower_[iSet];
          upper_[numberSets_] = upper_[iSet];
        }
      }
    }
  }
  delete[] array;
  numberSets_++; // adjust
  // Find type of gub
  firstGub_ = numberColumns + 1;
  lastGub_ = -1;
  for (i = 0; i < numberColumns; i++) {
    if (backward_[i] >= 0) {
      firstGub_ = std::min(firstGub_, i);
      lastGub_ = std::max(lastGub_, i);
    }
  }
  if (lastGub_ > 0)
    lastGub_++;
  gubType_ = 0;
  for (i = firstGub_; i < lastGub_; i++) {
    if (backward_[i] < 0) {
      gubType_ = 1;
      break;
    }
  }

  // Make sure key is feasible if only key in set
}
ClpGubMatrix::ClpGubMatrix(
  const CoinPackedMatrix &rhs,
  int numberRows, const int *whichRows,
  int numberColumns, const int *whichColumns)
  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
  , sumDualInfeasibilities_(0.0)
  , sumPrimalInfeasibilities_(0.0)
  , sumOfRelaxedDualInfeasibilities_(0.0)
  , sumOfRelaxedPrimalInfeasibilities_(0.0)
  , start_(NULL)
  , end_(NULL)
  , lower_(NULL)
  , upper_(NULL)
  , backward_(NULL)
  , backToPivotRow_(NULL)
  , changeCost_(NULL)
  , keyVariable_(NULL)
  , next_(NULL)
  , toIndex_(NULL)
  , fromIndex_(NULL)
  , numberDualInfeasibilities_(0)
  , numberPrimalInfeasibilities_(0)
  , numberSets_(0)
  , saveNumber_(0)
  , possiblePivotKey_(0)
  , gubSlackIn_(-1)
  , firstGub_(0)
  , lastGub_(0)
  , gubType_(0)
{
  setType(16);
}
/* Return <code>x * A + y</code> in <code>z</code>.
	Squashes small elements and knows about ClpSimplex */
void ClpGubMatrix::transposeTimes(const ClpSimplex *model, double scalar,
  const CoinIndexedVector *rowArray,
  CoinIndexedVector *y,
  CoinIndexedVector *columnArray) const
{
  columnArray->clear();
  double *pi = rowArray->denseVector();
  int numberNonZero = 0;
  int *index = columnArray->getIndices();
  double *array = columnArray->denseVector();
  int numberInRowArray = rowArray->getNumElements();
  // maybe I need one in OsiSimplex
  double zeroTolerance = model->zeroTolerance();
  int numberRows = model->numberRows();
  ClpPackedMatrix *rowCopy = dynamic_cast< ClpPackedMatrix * >(model->rowCopy());
  bool packed = rowArray->packedMode();
  double factor = 0.3;
  // We may not want to do by row if there may be cache problems
  int numberColumns = model->numberColumns();
  // It would be nice to find L2 cache size - for moment 512K
  // Be slightly optimistic
  if (numberColumns * sizeof(double) > 1000000) {
    if (numberRows * 10 < numberColumns)
      factor = 0.1;
    else if (numberRows * 4 < numberColumns)
      factor = 0.15;
    else if (numberRows * 2 < numberColumns)
      factor = 0.2;
    //if (model->numberIterations()%50==0)
    //printf("%d nonzero\n",numberInRowArray);
  }
  // reduce for gub
  factor *= 0.5;
  assert(!y->getNumElements());
  if (numberInRowArray > factor * numberRows || !rowCopy) {
    // do by column
    int iColumn;
    // get matrix data pointers
    const int *row = matrix_->getIndices();
    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
    const int *columnLength = matrix_->getVectorLengths();
    const double *elementByColumn = matrix_->getElements();
    const double *rowScale = model->rowScale();
    int numberColumns = model->numberColumns();
    int iSet = -1;
    double djMod = 0.0;
    if (packed) {
      // need to expand pi into y
      assert(y->capacity() >= numberRows);
      double *piOld = pi;
      pi = y->denseVector();
      const int *whichRow = rowArray->getIndices();
      int i;
      if (!rowScale) {
        // modify pi so can collapse to one loop
        for (i = 0; i < numberInRowArray; i++) {
          int iRow = whichRow[i];
          pi[iRow] = scalar * piOld[i];
        }
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          if (backward_[iColumn] != iSet) {
            // get pi on gub row
            iSet = backward_[iColumn];
            if (iSet >= 0) {
              int iBasic = keyVariable_[iSet];
              if (iBasic < numberColumns) {
                // get dj without
                assert(model->getStatus(iBasic) == ClpSimplex::basic);
                djMod = 0.0;
                for (CoinBigIndex j = columnStart[iBasic];
                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
                  int jRow = row[j];
                  djMod -= pi[jRow] * elementByColumn[j];
                }
              } else {
                djMod = 0.0;
              }
            } else {
              djMod = 0.0;
            }
          }
          double value = -djMod;
          CoinBigIndex j;
          for (j = columnStart[iColumn];
               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
            int iRow = row[j];
            value += pi[iRow] * elementByColumn[j];
          }
          if (fabs(value) > zeroTolerance) {
            array[numberNonZero] = value;
            index[numberNonZero++] = iColumn;
          }
        }
      } else {
        // scaled
        // modify pi so can collapse to one loop
        for (i = 0; i < numberInRowArray; i++) {
          int iRow = whichRow[i];
          pi[iRow] = scalar * piOld[i] * rowScale[iRow];
        }
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          if (backward_[iColumn] != iSet) {
            // get pi on gub row
            iSet = backward_[iColumn];
            if (iSet >= 0) {
              int iBasic = keyVariable_[iSet];
              if (iBasic < numberColumns) {
                // get dj without
                assert(model->getStatus(iBasic) == ClpSimplex::basic);
                djMod = 0.0;
                // scaled
                for (CoinBigIndex j = columnStart[iBasic];
                     j < columnStart[iBasic] + columnLength[iBasic]; j++) {
                  int jRow = row[j];
                  djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow];
                }
              } else {
                djMod = 0.0;
              }
            } else {
              djMod = 0.0;
            }
          }
          double value = -djMod;
          CoinBigIndex j;
          const double *columnScale = model->columnScale();
          for (j = columnStart[iColumn];
               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
            int iRow = row[j];
            value += pi[iRow] * elementByColumn[j];
          }
          value *= columnScale[iColumn];
          if (fabs(value) > zeroTolerance) {
            array[numberNonZero] = value;
            index[numberNonZero++] = iColumn;
          }
        }
      }
      // zero out
      for (i = 0; i < numberInRowArray; i++) {
        int iRow = whichRow[i];
        pi[iRow] = 0.0;
      }
    } else {
      // code later
      assert(packed);
      if (!rowScale) {
        if (scalar == -1.0) {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j];
            }
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = -value;
            }
          }
        } else if (scalar == 1.0) {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j];
            }
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = value;
            }
          }
        } else {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j];
            }
            value *= scalar;
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = value;
            }
          }
        }
      } else {
        // scaled
        if (scalar == -1.0) {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            const double *columnScale = model->columnScale();
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
            }
            value *= columnScale[iColumn];
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = -value;
            }
          }
        } else if (scalar == 1.0) {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            const double *columnScale = model->columnScale();
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
            }
            value *= columnScale[iColumn];
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = value;
            }
          }
        } else {
          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            double value = 0.0;
            CoinBigIndex j;
            const double *columnScale = model->columnScale();
            for (j = columnStart[iColumn];
                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
              int iRow = row[j];
              value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
            }
            value *= scalar * columnScale[iColumn];
            if (fabs(value) > zeroTolerance) {
              index[numberNonZero++] = iColumn;
              array[iColumn] = value;
            }
          }
        }
      }
    }
    columnArray->setNumElements(numberNonZero);
    y->setNumElements(0);
  } else {
    // do by row
    transposeTimesByRow(model, scalar, rowArray, y, columnArray);
  }
  if (packed)
    columnArray->setPackedMode(true);
  if (0) {
    columnArray->checkClean();
    int numberNonZero = columnArray->getNumElements();
    ;
    int *index = columnArray->getIndices();
    double *array = columnArray->denseVector();
    int i;
    for (i = 0; i < numberNonZero; i++) {
      int j = index[i];
      double value;
      if (packed)
        value = array[i];
      else
        value = array[j];
      printf("Ti %d %d %g\n", i, j, value);
    }
  }
}
/* Return <code>x * A + y</code> in <code>z</code>.
	Squashes small elements and knows about ClpSimplex */
void ClpGubMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
  const CoinIndexedVector *rowArray,
  CoinIndexedVector *y,
  CoinIndexedVector *columnArray) const
{
  // Do packed part
  ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
  if (numberSets_) {
    /* what we need to do is do by row as normal but get list of sets touched
             and then update those ones */
    abort();
  }
}
/* Return <code>x *A in <code>z</code> but
   just for indices in y. */
void ClpGubMatrix::subsetTransposeTimes(const ClpSimplex *model,
  const CoinIndexedVector *rowArray,
  const CoinIndexedVector *y,
  CoinIndexedVector *columnArray) const
{
  columnArray->clear();
  double *pi = rowArray->denseVector();
  double *array = columnArray->denseVector();
  int jColumn;
  // get matrix data pointers
  const int *row = matrix_->getIndices();
  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
  const int *columnLength = matrix_->getVectorLengths();
  const double *elementByColumn = matrix_->getElements();
  const double *rowScale = model->rowScale();
  int numberToDo = y->getNumElements();
  const int *which = y->getIndices();
  assert(!rowArray->packedMode());
  columnArray->setPacked();
  int numberTouched = 0;
  if (!rowScale) {
    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j;
      for (j = columnStart[iColumn];
           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
        int iRow = row[j];
        value += pi[iRow] * elementByColumn[j];
      }
      array[jColumn] = value;
      if (value) {
        int iSet = backward_[iColumn];
        if (iSet >= 0) {
          int iBasic = keyVariable_[iSet];
          if (iBasic == iColumn) {
            toIndex_[iSet] = jColumn;
            fromIndex_[numberTouched++] = iSet;
          }
        }
      }
    }
  } else {
    // scaled
    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j;
      const double *columnScale = model->columnScale();
      for (j = columnStart[iColumn];
           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
        int iRow = row[j];
        value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
      }
      value *= columnScale[iColumn];
      array[jColumn] = value;
      if (value) {
        int iSet = backward_[iColumn];
        if (iSet >= 0) {
          int iBasic = keyVariable_[iSet];
          if (iBasic == iColumn) {
            toIndex_[iSet] = jColumn;
            fromIndex_[numberTouched++] = iSet;
          }
        }
      }
    }
  }
  // adjust djs
  for (jColumn = 0; jColumn < numberToDo; jColumn++) {
    int iColumn = which[jColumn];
    int iSet = backward_[iColumn];
    if (iSet >= 0) {
      int kColumn = toIndex_[iSet];
      if (kColumn >= 0)
        array[jColumn] -= array[kColumn];
    }
  }
  // and clear basic
  for (int j = 0; j < numberTouched; j++) {
    int iSet = fromIndex_[j];
    int kColumn = toIndex_[iSet];
    toIndex_[iSet] = -1;
    array[kColumn] = 0.0;
  }
}
/// returns number of elements in column part of basis,
int ClpGubMatrix::countBasis(const int *whichColumn,
  int &numberColumnBasic)
{
  int i;
  int numberColumns = getNumCols();
  const int *columnLength = matrix_->getVectorLengths();
  int numberRows = getNumRows();
  int numberBasic = 0;
  CoinBigIndex numberElements = 0;
  int lastSet = -1;
  int key = -1;
  int keyLength = -1;
  double *work = new double[numberRows];
  CoinZeroN(work, numberRows);
  char *mark = new char[numberRows];
  CoinZeroN(mark, numberRows);
  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
  const int *row = matrix_->getIndices();
  const double *elementByColumn = matrix_->getElements();
  //ClpGubDynamicMatrix* gubx =
  //dynamic_cast< ClpGubDynamicMatrix*>(this);
  //int * id = gubx->id();
  // just count
  for (i = 0; i < numberColumnBasic; i++) {
    int iColumn = whichColumn[i];
    int iSet = backward_[iColumn];
    int length = columnLength[iColumn];
    if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
      numberElements += length;
      numberBasic++;
      //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
    } else {
      // in gub set
      if (iColumn != keyVariable_[iSet]) {
        numberBasic++;
        CoinBigIndex j;
        // not key
        if (lastSet < iSet) {
          // erase work
          if (key >= 0) {
            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
              work[row[j]] = 0.0;
          }
          key = keyVariable_[iSet];
          lastSet = iSet;
          keyLength = columnLength[key];
          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
            work[row[j]] = elementByColumn[j];
        }
        int extra = keyLength;
        for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
          int iRow = row[j];
          double keyValue = work[iRow];
          double value = elementByColumn[j];
          if (!keyValue) {
            if (fabs(value) > 1.0e-20)
              extra++;
          } else {
            value -= keyValue;
            if (fabs(value) <= 1.0e-20)
              extra--;
          }
        }
        numberElements += extra;
        //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
      }
    }
  }
  delete[] work;
  delete[] mark;
  // update number of column basic
  numberColumnBasic = numberBasic;
  if (numberElements > COIN_INT_MAX) {
    printf("Factorization too large\n");
    abort();
  }
  return static_cast< int >(numberElements);
}
void ClpGubMatrix::fillBasis(ClpSimplex *model,
  const int *whichColumn,
  int &numberColumnBasic,
  int *indexRowU, int *start,
  int *rowCount, int *columnCount,
  CoinFactorizationDouble *elementU)
{
  int i;
  int numberColumns = getNumCols();
  const int *columnLength = matrix_->getVectorLengths();
  int numberRows = getNumRows();
  assert(next_ || !elementU);
  CoinBigIndex numberElements = start[0];
  int lastSet = -1;
  int key = -1;
  int keyLength = -1;
  double *work = new double[numberRows];
  CoinZeroN(work, numberRows);
  char *mark = new char[numberRows];
  CoinZeroN(mark, numberRows);
  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
  const int *row = matrix_->getIndices();
  const double *elementByColumn = matrix_->getElements();
  const double *rowScale = model->rowScale();
  int numberBasic = 0;
  if (0) {
    printf("%d basiccolumns\n", numberColumnBasic);
    int i;
    for (i = 0; i < numberSets_; i++) {
      int k = keyVariable_[i];
      if (k < numberColumns) {
        printf("key %d on set %d, %d elements\n", k, i, columnStart[k + 1] - columnStart[k]);
        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
          printf("row %d el %g\n", row[j], elementByColumn[j]);
      } else {
        printf("slack key on set %d\n", i);
      }
    }
  }
  // fill
  if (!rowScale) {
    // no scaling
    for (i = 0; i < numberColumnBasic; i++) {
      int iColumn = whichColumn[i];
      int iSet = backward_[iColumn];
      int length = columnLength[iColumn];
      if (0) {
        int k = iColumn;
        printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k + 1] - columnStart[k]);
        for (int j = columnStart[k]; j < columnStart[k + 1]; j++)
          printf("row %d el %g\n", row[j], elementByColumn[j]);
      }
      CoinBigIndex j;
      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
          double value = elementByColumn[j];
          if (fabs(value) > 1.0e-20) {
            int iRow = row[j];
            indexRowU[numberElements] = iRow;
            rowCount[iRow]++;
            elementU[numberElements++] = value;
          }
        }
        // end of column
        columnCount[numberBasic] = numberElements - start[numberBasic];
        numberBasic++;
        start[numberBasic] = numberElements;
      } else {
        // in gub set
        if (iColumn != keyVariable_[iSet]) {
          // not key
          if (lastSet != iSet) {
            // erase work
            if (key >= 0) {
              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
                int iRow = row[j];
                work[iRow] = 0.0;
                mark[iRow] = 0;
              }
            }
            key = keyVariable_[iSet];
            lastSet = iSet;
            keyLength = columnLength[key];
            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
              int iRow = row[j];
              work[iRow] = elementByColumn[j];
              mark[iRow] = 1;
            }
          }
          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
            int iRow = row[j];
            double value = elementByColumn[j];
            if (mark[iRow]) {
              mark[iRow] = 0;
              double keyValue = work[iRow];
              value -= keyValue;
            }
            if (fabs(value) > 1.0e-20) {
              indexRowU[numberElements] = iRow;
              rowCount[iRow]++;
              elementU[numberElements++] = value;
            }
          }
          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
            int iRow = row[j];
            if (mark[iRow]) {
              double value = -work[iRow];
              if (fabs(value) > 1.0e-20) {
                indexRowU[numberElements] = iRow;
                rowCount[iRow]++;
                elementU[numberElements++] = value;
              }
            } else {
              // just put back mark
              mark[iRow] = 1;
            }
          }
          // end of column
          columnCount[numberBasic] = numberElements - start[numberBasic];
          numberBasic++;
          start[numberBasic] = numberElements;
        }
      }
    }
  } else {
    // scaling
    const double *columnScale = model->columnScale();
    for (i = 0; i < numberColumnBasic; i++) {
      int iColumn = whichColumn[i];
      int iSet = backward_[iColumn];
      int length = columnLength[iColumn];
      CoinBigIndex j;
      if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
        double scale = columnScale[iColumn];
        for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
          int iRow = row[j];
          double value = elementByColumn[j] * scale * rowScale[iRow];
          if (fabs(value) > 1.0e-20) {
            indexRowU[numberElements] = iRow;
            rowCount[iRow]++;
            elementU[numberElements++] = value;
          }
        }
        // end of column
        columnCount[numberBasic] = numberElements - start[numberBasic];
        numberBasic++;
        start[numberBasic] = numberElements;
      } else {
        // in gub set
        if (iColumn != keyVariable_[iSet]) {
          double scale = columnScale[iColumn];
          // not key
          if (lastSet < iSet) {
            // erase work
            if (key >= 0) {
              for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
                int iRow = row[j];
                work[iRow] = 0.0;
                mark[iRow] = 0;
              }
            }
            key = keyVariable_[iSet];
            lastSet = iSet;
            keyLength = columnLength[key];
            double scale = columnScale[key];
            for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
              int iRow = row[j];
              work[iRow] = elementByColumn[j] * scale * rowScale[iRow];
              mark[iRow] = 1;
            }
          }
          for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
            int iRow = row[j];
            double value = elementByColumn[j] * scale * rowScale[iRow];
            if (mark[iRow]) {
              mark[iRow] = 0;
              double keyValue = work[iRow];
              value -= keyValue;
            }
            if (fabs(value) > 1.0e-20) {
              indexRowU[numberElements] = iRow;
              rowCount[iRow]++;
              elementU[numberElements++] = value;
            }
          }
          for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
            int iRow = row[j];
            if (mark[iRow]) {
              double value = -work[iRow];
              if (fabs(value) > 1.0e-20) {
                indexRowU[numberElements] = iRow;
                rowCount[iRow]++;
                elementU[numberElements++] = value;
              }
            } else {
              // just put back mark
              mark[iRow] = 1;
            }
          }
          // end of column
          columnCount[numberBasic] = numberElements - start[numberBasic];
          numberBasic++;
          start[numberBasic] = numberElements;
        }
      }
    }
  }
  delete[] work;
  delete[] mark;
  // update number of column basic
  numberColumnBasic = numberBasic;
}
/* Unpacks a column into an CoinIndexedvector
 */
void ClpGubMatrix::unpack(const ClpSimplex *model, CoinIndexedVector *rowArray,
  int iColumn) const
{
  assert(iColumn < model->numberColumns());
  // Do packed part
  ClpPackedMatrix::unpack(model, rowArray, iColumn);
  int iSet = backward_[iColumn];
  if (iSet >= 0) {
    int iBasic = keyVariable_[iSet];
    if (iBasic < model->numberColumns()) {
      add(model, rowArray, iBasic, -1.0);
    }
  }
}
/* Unpacks a column into a CoinIndexedVector
** in packed format
Note that model is NOT const.  Bounds and objective could
be modified if doing column generation (just for this variable) */
void ClpGubMatrix::unpackPacked(ClpSimplex *model,
  CoinIndexedVector *rowArray,
  int iColumn) const
{
  int numberColumns = model->numberColumns();
  if (iColumn < numberColumns) {
    // Do packed part
    ClpPackedMatrix::unpackPacked(model, rowArray, iColumn);
    int iSet = backward_[iColumn];
    if (iSet >= 0) {
      // columns are in order
      int iBasic = keyVariable_[iSet];
      if (iBasic < numberColumns) {
        int number = rowArray->getNumElements();
        const double *rowScale = model->rowScale();
        const int *row = matrix_->getIndices();
        const CoinBigIndex *columnStart = matrix_->getVectorStarts();
        const int *columnLength = matrix_->getVectorLengths();
        const double *elementByColumn = matrix_->getElements();
        double *array = rowArray->denseVector();
        int *index = rowArray->getIndices();
        CoinBigIndex i;
        int numberOld = number;
        int lastIndex = 0;
        int next = index[lastIndex];
        if (!rowScale) {
          for (i = columnStart[iBasic];
               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
            int iRow = row[i];
            while (iRow > next) {
              lastIndex++;
              if (lastIndex == numberOld)
                next = matrix_->getNumRows();
              else
                next = index[lastIndex];
            }
            if (iRow < next) {
              array[number] = -elementByColumn[i];
              index[number++] = iRow;
            } else {
              assert(iRow == next);
              array[lastIndex] -= elementByColumn[i];
              if (!array[lastIndex])
                array[lastIndex] = 1.0e-100;
            }
          }
        } else {
          // apply scaling
          double scale = model->columnScale()[iBasic];
          for (i = columnStart[iBasic];
               i < columnStart[iBasic] + columnLength[iBasic]; i++) {
            int iRow = row[i];
            while (iRow > next) {
              lastIndex++;
              if (lastIndex == numberOld)
                next = matrix_->getNumRows();
              else
                next = index[lastIndex];
            }
            if (iRow < next) {
              array[number] = -elementByColumn[i] * scale * rowScale[iRow];
              index[number++] = iRow;
            } else {
              assert(iRow == next);
              array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow];
              if (!array[lastIndex])
                array[lastIndex] = 1.0e-100;
            }
          }
        }
        rowArray->setNumElements(number);
      }
    }
  } else {
    // key slack entering
    int iBasic = keyVariable_[gubSlackIn_];
    assert(iBasic < numberColumns);
    int number = 0;
    const double *rowScale = model->rowScale();
    const int *row = matrix_->getIndices();
    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
    const int *columnLength = matrix_->getVectorLengths();
    const double *elementByColumn = matrix_->getElements();
    double *array = rowArray->denseVector();
    int *index = rowArray->getIndices();
    CoinBigIndex i;
    if (!rowScale) {
      for (i = columnStart[iBasic];
           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
        int iRow = row[i];
        array[number] = elementByColumn[i];
        index[number++] = iRow;
      }
    } else {
      // apply scaling
      double scale = model->columnScale()[iBasic];
      for (i = columnStart[iBasic];
           i < columnStart[iBasic] + columnLength[iBasic]; i++) {
        int iRow = row[i];
        array[number] = elementByColumn[i] * scale * rowScale[iRow];
        index[number++] = iRow;
      }
    }
    rowArray->setNumElements(number);
    rowArray->setPacked();
  }
}
/* Adds multiple of a column into an CoinIndexedvector
      You can use quickAdd to add to vector */
void ClpGubMatrix::add(const ClpSimplex *model, CoinIndexedVector *rowArray,
  int iColumn, double multiplier) const
{
  assert(iColumn < model->numberColumns());
  // Do packed part
  ClpPackedMatrix::add(model, rowArray, iColumn, multiplier);
  int iSet = backward_[iColumn];
  if (iSet >= 0 && iColumn != keyVariable_[iSet]) {
    ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier);
  }
}
/* Adds multiple of a column into an array */
void ClpGubMatrix::add(const ClpSimplex *model, double *array,
  int iColumn, double multiplier) const
{
  assert(iColumn < model->numberColumns());
  // Do packed part
  ClpPackedMatrix::add(model, array, iColumn, multiplier);
  if (iColumn < model->numberColumns()) {
    int iSet = backward_[iColumn];
    if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) {
      ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier);
    }
  }
}
// Partial pricing
void ClpGubMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
  int &bestSequence, int &numberWanted)
{
  numberWanted = currentWanted_;
  if (numberSets_) {
    // Do packed part before gub
    int numberColumns = matrix_->getNumCols();
    double ratio = static_cast< double >(firstGub_) / static_cast< double >(numberColumns);
    ClpPackedMatrix::partialPricing(model, startFraction * ratio,
      endFraction * ratio, bestSequence, numberWanted);
    if (numberWanted || minimumGoodReducedCosts_ < -1) {
      // do gub
      const double *element = matrix_->getElements();
      const int *row = matrix_->getIndices();
      const CoinBigIndex *startColumn = matrix_->getVectorStarts();
      const int *length = matrix_->getVectorLengths();
      const double *rowScale = model->rowScale();
      const double *columnScale = model->columnScale();
      int iSequence;
      CoinBigIndex j;
      double tolerance = model->currentDualTolerance();
      double *reducedCost = model->djRegion();
      const double *duals = model->dualRowSolution();
      const double *cost = model->costRegion();
      double bestDj;
      int numberColumns = model->numberColumns();
      int numberRows = model->numberRows();
      if (bestSequence >= 0)
        bestDj = fabs(this->reducedCost(model, bestSequence));
      else
        bestDj = tolerance;
      int sequenceOut = model->sequenceOut();
      int saveSequence = bestSequence;
      int startG = firstGub_ + static_cast< int >(startFraction * (lastGub_ - firstGub_));
      int endG = firstGub_ + static_cast< int >(endFraction * (lastGub_ - firstGub_));
      endG = std::min(lastGub_, endG + 1);
      // If nothing found yet can go all the way to end
      int endAll = endG;
      if (bestSequence < 0 && !startG)
        endAll = lastGub_;
      int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
      int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_;
      int nSets = 0;
      int iSet = -1;
      double djMod = 0.0;
      double infeasibilityCost = model->infeasibilityCost();
      if (rowScale) {
        double bestDjMod = 0.0;
        // scaled
        for (iSequence = startG; iSequence < endAll; iSequence++) {
          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
            // give up
            numberWanted = 0;
            break;
          } else if (iSequence == endG && bestSequence >= 0) {
            break;
          }
          if (backward_[iSequence] != iSet) {
            // get pi on gub row
            iSet = backward_[iSequence];
            if (iSet >= 0) {
              nSets++;
              int iBasic = keyVariable_[iSet];
              if (iBasic >= numberColumns) {
                djMod = -weight(iSet) * infeasibilityCost;
              } else {
                // get dj without
                assert(model->getStatus(iBasic) == ClpSimplex::basic);
                djMod = 0.0;
                // scaled
                for (j = startColumn[iBasic];
                     j < startColumn[iBasic] + length[iBasic]; j++) {
                  int jRow = row[j];
                  djMod -= duals[jRow] * element[j] * rowScale[jRow];
                }
                // allow for scaling
                djMod += cost[iBasic] / columnScale[iBasic];
                // See if gub slack possible - dj is djMod
                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
                  double value = -djMod;
                  if (value > tolerance) {
                    numberWanted--;
                    if (value > bestDj) {
                      // check flagged variable and correct dj
                      if (!flagged(iSet)) {
                        bestDj = value;
                        bestSequence = numberRows + numberColumns + iSet;
                        bestDjMod = djMod;
                      } else {
                        // just to make sure we don't exit before got something
                        numberWanted++;
                        abort();
                      }
                    }
                  }
                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
                  double value = djMod;
                  if (value > tolerance) {
                    numberWanted--;
                    if (value > bestDj) {
                      // check flagged variable and correct dj
                      if (!flagged(iSet)) {
                        bestDj = value;
                        bestSequence = numberRows + numberColumns + iSet;
                        bestDjMod = djMod;
                      } else {
                        // just to make sure we don't exit before got something
                        numberWanted++;
                        abort();
                      }
                    }
                  }
                }
              }
            } else {
              // not in set
              djMod = 0.0;
            }
          }
          if (iSequence != sequenceOut) {
            double value;
            ClpSimplex::Status status = model->getStatus(iSequence);

            switch (status) {

            case ClpSimplex::basic:
            case ClpSimplex::isFixed:
              break;
            case ClpSimplex::isFree:
            case ClpSimplex::superBasic:
              value = -djMod;
              // scaled
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j] * rowScale[jRow];
              }
              value = fabs(cost[iSequence] + value * columnScale[iSequence]);
              if (value > FREE_ACCEPT * tolerance) {
                numberWanted--;
                // we are going to bias towards free (but only if reasonable)
                value *= FREE_BIAS;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            case ClpSimplex::atUpperBound:
              value = -djMod;
              // scaled
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j] * rowScale[jRow];
              }
              value = cost[iSequence] + value * columnScale[iSequence];
              if (value > tolerance) {
                numberWanted--;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            case ClpSimplex::atLowerBound:
              value = -djMod;
              // scaled
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j] * rowScale[jRow];
              }
              value = -(cost[iSequence] + value * columnScale[iSequence]);
              if (value > tolerance) {
                numberWanted--;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            }
          }
          if (!numberWanted)
            break;
        }
        if (bestSequence != saveSequence) {
          if (bestSequence < numberRows + numberColumns) {
            // recompute dj
            double value = bestDjMod;
            // scaled
            for (j = startColumn[bestSequence];
                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
              int jRow = row[j];
              value -= duals[jRow] * element[j] * rowScale[jRow];
            }
            reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
            gubSlackIn_ = -1;
          } else {
            // slack - make last column
            gubSlackIn_ = bestSequence - numberRows - numberColumns;
            bestSequence = numberColumns + 2 * numberRows;
            reducedCost[bestSequence] = bestDjMod;
            model->setStatus(bestSequence, getStatus(gubSlackIn_));
            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
            else
              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
            model->costRegion()[bestSequence] = 0.0;
          }
          savedBestSequence_ = bestSequence;
          savedBestDj_ = reducedCost[savedBestSequence_];
        }
      } else {
        double bestDjMod = 0.0;
        //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
        //     startG,endG,numberWanted);
        for (iSequence = startG; iSequence < endG; iSequence++) {
          if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
            // give up
            numberWanted = 0;
            break;
          } else if (iSequence == endG && bestSequence >= 0) {
            break;
          }
          if (backward_[iSequence] != iSet) {
            // get pi on gub row
            iSet = backward_[iSequence];
            if (iSet >= 0) {
              nSets++;
              int iBasic = keyVariable_[iSet];
              if (iBasic >= numberColumns) {
                djMod = -weight(iSet) * infeasibilityCost;
              } else {
                // get dj without
                assert(model->getStatus(iBasic) == ClpSimplex::basic);
                djMod = 0.0;

                for (j = startColumn[iBasic];
                     j < startColumn[iBasic] + length[iBasic]; j++) {
                  int jRow = row[j];
                  djMod -= duals[jRow] * element[j];
                }
                djMod += cost[iBasic];
                // See if gub slack possible - dj is djMod
                if (getStatus(iSet) == ClpSimplex::atLowerBound) {
                  double value = -djMod;
                  if (value > tolerance) {
                    numberWanted--;
                    if (value > bestDj) {
                      // check flagged variable and correct dj
                      if (!flagged(iSet)) {
                        bestDj = value;
                        bestSequence = numberRows + numberColumns + iSet;
                        bestDjMod = djMod;
                      } else {
                        // just to make sure we don't exit before got something
                        numberWanted++;
                        abort();
                      }
                    }
                  }
                } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
                  double value = djMod;
                  if (value > tolerance) {
                    numberWanted--;
                    if (value > bestDj) {
                      // check flagged variable and correct dj
                      if (!flagged(iSet)) {
                        bestDj = value;
                        bestSequence = numberRows + numberColumns + iSet;
                        bestDjMod = djMod;
                      } else {
                        // just to make sure we don't exit before got something
                        numberWanted++;
                        abort();
                      }
                    }
                  }
                }
              }
            } else {
              // not in set
              djMod = 0.0;
            }
          }
          if (iSequence != sequenceOut) {
            double value;
            ClpSimplex::Status status = model->getStatus(iSequence);

            switch (status) {

            case ClpSimplex::basic:
            case ClpSimplex::isFixed:
              break;
            case ClpSimplex::isFree:
            case ClpSimplex::superBasic:
              value = cost[iSequence] - djMod;
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j];
              }
              value = fabs(value);
              if (value > FREE_ACCEPT * tolerance) {
                numberWanted--;
                // we are going to bias towards free (but only if reasonable)
                value *= FREE_BIAS;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            case ClpSimplex::atUpperBound:
              value = cost[iSequence] - djMod;
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j];
              }
              if (value > tolerance) {
                numberWanted--;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            case ClpSimplex::atLowerBound:
              value = cost[iSequence] - djMod;
              for (j = startColumn[iSequence];
                   j < startColumn[iSequence] + length[iSequence]; j++) {
                int jRow = row[j];
                value -= duals[jRow] * element[j];
              }
              value = -value;
              if (value > tolerance) {
                numberWanted--;
                if (value > bestDj) {
                  // check flagged variable and correct dj
                  if (!model->flagged(iSequence)) {
                    bestDj = value;
                    bestSequence = iSequence;
                    bestDjMod = djMod;
                  } else {
                    // just to make sure we don't exit before got something
                    numberWanted++;
                  }
                }
              }
              break;
            }
          }
          if (!numberWanted)
            break;
        }
        if (bestSequence != saveSequence) {
          if (bestSequence < numberRows + numberColumns) {
            // recompute dj
            double value = cost[bestSequence] - bestDjMod;
            for (j = startColumn[bestSequence];
                 j < startColumn[bestSequence] + length[bestSequence]; j++) {
              int jRow = row[j];
              value -= duals[jRow] * element[j];
            }
            //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
            reducedCost[bestSequence] = value;
            gubSlackIn_ = -1;
          } else {
            // slack - make last column
            gubSlackIn_ = bestSequence - numberRows - numberColumns;
            bestSequence = numberColumns + 2 * numberRows;
            reducedCost[bestSequence] = bestDjMod;
            //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
            model->setStatus(bestSequence, getStatus(gubSlackIn_));
            if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
              model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
            else
              model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
            model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
            model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
            model->costRegion()[bestSequence] = 0.0;
          }
        }
      }
      // See if may be finished
      if (startG == firstGub_ && bestSequence < 0)
        infeasibilityWeight_ = model_->infeasibilityCost();
      else if (bestSequence >= 0)
        infeasibilityWeight_ = -1.0;
    }
    if (numberWanted) {
      // Do packed part after gub
      double offset = static_cast< double >(lastGub_) / static_cast< double >(numberColumns);
      double ratio = static_cast< double >(numberColumns) / static_cast< double >(numberColumns) - offset;
      double start2 = offset + ratio * startFraction;
      double end2 = std::min(1.0, offset + ratio * endFraction + 1.0e-6);
      ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted);
    }
  } else {
    // no gub
    ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
  }
  if (bestSequence >= 0)
    infeasibilityWeight_ = -1.0; // not optimal
  currentWanted_ = numberWanted;
}
/* expands an updated column to allow for extra rows which the main
   solver does not know about and returns number added.
*/
int ClpGubMatrix::extendUpdated(ClpSimplex *model, CoinIndexedVector *update, int mode)
{
  // I think we only need to bother about sets with two in basis or incoming set
  int number = update->getNumElements();
  double *array = update->denseVector();
  int *index = update->getIndices();
  int i;
  assert(!number || update->packedMode());
  int *pivotVariable = model->pivotVariable();
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  int numberTotal = numberRows + numberColumns;
  int sequenceIn = model->sequenceIn();
  int returnCode = 0;
  int iSetIn;
  if (sequenceIn < numberColumns) {
    iSetIn = backward_[sequenceIn];
    gubSlackIn_ = -1; // in case set
  } else if (sequenceIn < numberRows + numberColumns) {
    iSetIn = -1;
    gubSlackIn_ = -1; // in case set
  } else {
    iSetIn = gubSlackIn_;
  }
  double *lower = model->lowerRegion();
  double *upper = model->upperRegion();
  double *cost = model->costRegion();
  double *solution = model->solutionRegion();
  int number2 = number;
  if (!mode) {
    double primalTolerance = model->primalTolerance();
    double infeasibilityCost = model->infeasibilityCost();
    // extend
    saveNumber_ = number;
    for (i = 0; i < number; i++) {
      int iRow = index[i];
      int iPivot = pivotVariable[iRow];
      if (iPivot < numberColumns) {
        int iSet = backward_[iPivot];
        if (iSet >= 0) {
          // two (or more) in set
          int iIndex = toIndex_[iSet];
          double otherValue = array[i];
          double value;
          if (iIndex < 0) {
            toIndex_[iSet] = number2;
            int iNew = number2 - number;
            fromIndex_[number2 - number] = iSet;
            iIndex = number2;
            index[number2] = numberRows + iNew;
            // do key stuff
            int iKey = keyVariable_[iSet];
            if (iKey < numberColumns) {
              // Save current cost of key
              changeCost_[number2 - number] = cost[iKey];
              if (iSet != iSetIn)
                value = 0.0;
              else if (iSetIn != gubSlackIn_)
                value = 1.0;
              else
                value = -1.0;
              pivotVariable[numberRows + iNew] = iKey;
              // Do I need to recompute?
              double sol;
              assert(getStatus(iSet) != ClpSimplex::basic);
              if (getStatus(iSet) == ClpSimplex::atLowerBound)
                sol = lower_[iSet];
              else
                sol = upper_[iSet];
              if ((gubType_ & 8) != 0) {
                int iColumn = next_[iKey];
                // sum all non-key variables
                while (iColumn >= 0) {
                  sol -= solution[iColumn];
                  iColumn = next_[iColumn];
                }
              } else {
                int stop = -(iKey + 1);
                int iColumn = next_[iKey];
                // sum all non-key variables
                while (iColumn != stop) {
                  if (iColumn < 0)
                    iColumn = -iColumn - 1;
                  sol -= solution[iColumn];
                  iColumn = next_[iColumn];
                }
              }
              solution[iKey] = sol;
              if (model->algorithm() > 0)
                model->nonLinearCost()->setOne(iKey, sol);
              //assert (fabs(sol-solution[iKey])<1.0e-3);
            } else {
              // gub slack is basic
              // Save current cost of key
              changeCost_[number2 - number] = -weight(iSet) * infeasibilityCost;
              otherValue = -otherValue; //allow for - sign on slack
              if (iSet != iSetIn)
                value = 0.0;
              else
                value = -1.0;
              pivotVariable[numberRows + iNew] = iNew + numberTotal;
              model->djRegion()[iNew + numberTotal] = 0.0;
              double sol = 0.0;
              if ((gubType_ & 8) != 0) {
                int iColumn = next_[iKey];
                // sum all non-key variables
                while (iColumn >= 0) {
                  sol += solution[iColumn];
                  iColumn = next_[iColumn];
                }
              } else {
                int stop = -(iKey + 1);
                int iColumn = next_[iKey];
                // sum all non-key variables
                while (iColumn != stop) {
                  if (iColumn < 0)
                    iColumn = -iColumn - 1;
                  sol += solution[iColumn];
                  iColumn = next_[iColumn];
                }
              }
              solution[iNew + numberTotal] = sol;
              // and do cost in nonLinearCost
              if (model->algorithm() > 0)
                model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]);
              if (sol > upper_[iSet] + primalTolerance) {
                setAbove(iSet);
                lower[iNew + numberTotal] = upper_[iSet];
                upper[iNew + numberTotal] = COIN_DBL_MAX;
              } else if (sol < lower_[iSet] - primalTolerance) {
                setBelow(iSet);
                lower[iNew + numberTotal] = -COIN_DBL_MAX;
                upper[iNew + numberTotal] = lower_[iSet];
              } else {
                setFeasible(iSet);
                lower[iNew + numberTotal] = lower_[iSet];
                upper[iNew + numberTotal] = upper_[iSet];
              }
              cost[iNew + numberTotal] = weight(iSet) * infeasibilityCost;
            }
            number2++;
          } else {
            value = array[iIndex];
            int iKey = keyVariable_[iSet];
            if (iKey >= numberColumns)
              otherValue = -otherValue; //allow for - sign on slack
          }
          value -= otherValue;
          array[iIndex] = value;
        }
      }
    }
    if (iSetIn >= 0 && toIndex_[iSetIn] < 0) {
      // Do incoming
      update->setPacked(); // just in case no elements
      toIndex_[iSetIn] = number2;
      int iNew = number2 - number;
      fromIndex_[number2 - number] = iSetIn;
      // Save current cost of key
      double currentCost;
      int key = keyVariable_[iSetIn];
      if (key < numberColumns)
        currentCost = cost[key];
      else
        currentCost = -weight(iSetIn) * infeasibilityCost;
      changeCost_[number2 - number] = currentCost;
      index[number2] = numberRows + iNew;
      // do key stuff
      int iKey = keyVariable_[iSetIn];
      if (iKey < numberColumns) {
        if (gubSlackIn_ < 0)
          array[number2] = 1.0;
        else
          array[number2] = -1.0;
        pivotVariable[numberRows + iNew] = iKey;
        // Do I need to recompute?
        double sol;
        assert(getStatus(iSetIn) != ClpSimplex::basic);
        if (getStatus(iSetIn) == ClpSimplex::atLowerBound)
          sol = lower_[iSetIn];
        else
          sol = upper_[iSetIn];
        if ((gubType_ & 8) != 0) {
          int iColumn = next_[iKey];
          // sum all non-key variables
          while (iColumn >= 0) {
            sol -= solution[iColumn];
            iColumn = next_[iColumn];
          }
        } else {
          // bounds exist - sum over all except key
          int stop = -(iKey + 1);
          int iColumn = next_[iKey];
          // sum all non-key variables
          while (iColumn != stop) {
            if (iColumn < 0)
              iColumn = -iColumn - 1;
            sol -= solution[iColumn];
            iColumn = next_[iColumn];
          }
        }
        solution[iKey] = sol;
        if (model->algorithm() > 0)
          model->nonLinearCost()->setOne(iKey, sol);
        //assert (fabs(sol-solution[iKey])<1.0e-3);
      } else {
        // gub slack is basic
        array[number2] = -1.0;
        pivotVariable[numberRows + iNew] = iNew + numberTotal;
        model->djRegion()[iNew + numberTotal] = 0.0;
        double sol = 0.0;
        if ((gubType_ & 8) != 0) {
          int iColumn = next_[iKey];
          // sum all non-key variables
          while (iColumn >= 0) {
            sol += solution[iColumn];
            iColumn = next_[iColumn];
          }
        } else {
          // bounds exist - sum over all except key
          int stop = -(iKey + 1);
          int iColumn = next_[iKey];
          // sum all non-key variables
          while (iColumn != stop) {
            if (iColumn < 0)
              iColumn = -iColumn - 1;
            sol += solution[iColumn];
            iColumn = next_[iColumn];
          }
        }
        solution[iNew + numberTotal] = sol;
        // and do cost in nonLinearCost
        if (model->algorithm() > 0)
          model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]);
        if (sol > upper_[iSetIn] + primalTolerance) {
          setAbove(iSetIn);
          lower[iNew + numberTotal] = upper_[iSetIn];
          upper[iNew + numberTotal] = COIN_DBL_MAX;
        } else if (sol < lower_[iSetIn] - primalTolerance) {
          setBelow(iSetIn);
          lower[iNew + numberTotal] = -COIN_DBL_MAX;
          upper[iNew + numberTotal] = lower_[iSetIn];
        } else {
          setFeasible(iSetIn);
          lower[iNew + numberTotal] = lower_[iSetIn];
          upper[iNew + numberTotal] = upper_[iSetIn];
        }
        cost[iNew + numberTotal] = weight(iSetIn) * infeasibilityCost;
      }
      number2++;
    }
    // mark end
    fromIndex_[number2 - number] = -1;
    returnCode = number2 - number;
    // make sure lower_ upper_ adjusted
    synchronize(model, 9);
  } else {
    // take off?
    if (number > saveNumber_) {
      // clear
      double theta = model->theta();
      double *solution = model->solutionRegion();
      for (i = saveNumber_; i < number; i++) {
        int iRow = index[i];
        int iColumn = pivotVariable[iRow];
#ifdef CLP_DEBUG_PRINT
        printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
          iColumn, fromIndex_[i - saveNumber_], lower[iColumn], upper[iColumn], array[i],
          solution[iColumn], solution[iColumn] - model->theta() * array[i], model->theta());
#endif
        double value = array[i];
        array[i] = 0.0;
        int iSet = fromIndex_[i - saveNumber_];
        toIndex_[iSet] = -1;
        if (iSet == iSetIn && iColumn < numberColumns) {
          // update as may need value
          solution[iColumn] -= theta * value;
        }
      }
    }
#ifdef CLP_DEBUG
    for (i = 0; i < numberSets_; i++)
      assert(toIndex_[i] == -1);
#endif
    number2 = saveNumber_;
  }
  update->setNumElements(number2);
  return returnCode;
}
/*
     utility primal function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void ClpGubMatrix::primalExpanded(ClpSimplex *model, int mode)
{
  int numberColumns = model->numberColumns();
  switch (mode) {
    // If key variable then slot in gub rhs so will get correct contribution
  case 0: {
    int i;
    double *solution = model->solutionRegion();
    ClpSimplex::Status iStatus;
    for (i = 0; i < numberSets_; i++) {
      int iColumn = keyVariable_[i];
      if (iColumn < numberColumns) {
        // key is structural - where is slack
        iStatus = getStatus(i);
        assert(iStatus != ClpSimplex::basic);
        if (iStatus == ClpSimplex::atLowerBound)
          solution[iColumn] = lower_[i];
        else
          solution[iColumn] = upper_[i];
      }
    }
  } break;
  // Compute values of key variables
  case 1: {
    int i;
    double *solution = model->solutionRegion();
    //const int * columnLength = matrix_->getVectorLengths();
    //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    //const int * row = matrix_->getIndices();
    //const double * elementByColumn = matrix_->getElements();
    //int * pivotVariable = model->pivotVariable();
    sumPrimalInfeasibilities_ = 0.0;
    numberPrimalInfeasibilities_ = 0;
    double primalTolerance = model->primalTolerance();
    double relaxedTolerance = primalTolerance;
    // we can't really trust infeasibilities if there is primal error
    double error = std::min(1.0e-2, model->largestPrimalError());
    // allow tolerance at least slightly bigger than standard
    relaxedTolerance = relaxedTolerance + error;
    // but we will be using difference
    relaxedTolerance -= primalTolerance;
    sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
      int kColumn = keyVariable_[i];
      double value = 0.0;
      if ((gubType_ & 8) != 0) {
        int iColumn = next_[kColumn];
        // sum all non-key variables
        while (iColumn >= 0) {
          value += solution[iColumn];
          iColumn = next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(kColumn + 1);
        int iColumn = next_[kColumn];
        // sum all non-key variables
        while (iColumn != stop) {
          if (iColumn < 0)
            iColumn = -iColumn - 1;
          value += solution[iColumn];
          iColumn = next_[iColumn];
        }
      }
      if (kColumn < numberColumns) {
        // make sure key is basic - so will be skipped in values pass
        model->setStatus(kColumn, ClpSimplex::basic);
        // feasibility will be done later
        assert(getStatus(i) != ClpSimplex::basic);
        if (getStatus(i) == ClpSimplex::atUpperBound)
          solution[kColumn] = upper_[i] - value;
        else
          solution[kColumn] = lower_[i] - value;
        //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
      } else {
        // slack is key
        assert(getStatus(i) == ClpSimplex::basic);
        double infeasibility = 0.0;
        if (value > upper_[i] + primalTolerance) {
          infeasibility = value - upper_[i] - primalTolerance;
          setAbove(i);
        } else if (value < lower_[i] - primalTolerance) {
          infeasibility = lower_[i] - value - primalTolerance;
          setBelow(i);
        } else {
          setFeasible(i);
        }
        //printf("Value of key slack for set %d is %g\n",i,value);
        if (infeasibility > 0.0) {
          sumPrimalInfeasibilities_ += infeasibility;
          if (infeasibility > relaxedTolerance)
            sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
          numberPrimalInfeasibilities_++;
        }
      }
    }
  } break;
  // Report on infeasibilities of key variables
  case 2: {
    model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + sumPrimalInfeasibilities_);
    model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + numberPrimalInfeasibilities_);
    model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + sumOfRelaxedPrimalInfeasibilities_);
  } break;
  }
}
/*
     utility dual function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void ClpGubMatrix::dualExpanded(ClpSimplex *model,
  CoinIndexedVector *array,
  double * /*other*/, int mode)
{
  switch (mode) {
    // modify costs before transposeUpdate
  case 0: {
    int i;
    double *cost = model->costRegion();
    // not dual values yet
    //assert (!other);
    //double * work = array->denseVector();
    double infeasibilityCost = model->infeasibilityCost();
    int *pivotVariable = model->pivotVariable();
    int numberRows = model->numberRows();
    int numberColumns = model->numberColumns();
    for (i = 0; i < numberRows; i++) {
      int iPivot = pivotVariable[i];
      if (iPivot < numberColumns) {
        int iSet = backward_[iPivot];
        if (iSet >= 0) {
          int kColumn = keyVariable_[iSet];
          double costValue;
          if (kColumn < numberColumns) {
            // structural has cost
            costValue = cost[kColumn];
          } else {
            // slack is key
            assert(getStatus(iSet) == ClpSimplex::basic);
            // negative as -1.0 for slack
            costValue = -weight(iSet) * infeasibilityCost;
          }
          array->add(i, -costValue); // was work[i]-costValue
        }
      }
    }
  } break;
  // create duals for key variables (without check on dual infeasible)
  case 1: {
    // If key slack then dual 0.0 (if feasible)
    // dj for key is zero so that defines dual on set
    int i;
    double *dj = model->djRegion();
    int numberColumns = model->numberColumns();
    double infeasibilityCost = model->infeasibilityCost();
    for (i = 0; i < numberSets_; i++) {
      int kColumn = keyVariable_[i];
      if (kColumn < numberColumns) {
        // dj without set
        double value = dj[kColumn];
        // Now subtract out from all
        dj[kColumn] = 0.0;
        int iColumn = next_[kColumn];
        // modify all non-key variables
        while (iColumn >= 0) {
          dj[iColumn] -= value;
          iColumn = next_[iColumn];
        }
      } else {
        // slack key - may not be feasible
        assert(getStatus(i) == ClpSimplex::basic);
        // negative as -1.0 for slack
        double value = -weight(i) * infeasibilityCost;
        if (value) {
          int iColumn = next_[kColumn];
          // modify all non-key variables basic
          while (iColumn >= 0) {
            dj[iColumn] -= value;
            iColumn = next_[iColumn];
          }
        }
      }
    }
  } break;
  // as 1 but check slacks and compute djs
  case 2: {
    // If key slack then dual 0.0
    // If not then slack could be dual infeasible
    // dj for key is zero so that defines dual on set
    int i;
    // make sure fromIndex will not confuse pricing
    fromIndex_[0] = -1;
    possiblePivotKey_ = -1;
    // Create array
    int numberColumns = model->numberColumns();
    int *pivotVariable = model->pivotVariable();
    int numberRows = model->numberRows();
    for (i = 0; i < numberRows; i++) {
      int iPivot = pivotVariable[i];
      if (iPivot < numberColumns)
        backToPivotRow_[iPivot] = i;
    }
    if (noCheck_ >= 0) {
      if (infeasibilityWeight_ != model->infeasibilityCost()) {
        // don't bother checking
        sumDualInfeasibilities_ = 100.0;
        numberDualInfeasibilities_ = 1;
        sumOfRelaxedDualInfeasibilities_ = 100.0;
        return;
      }
    }
    double *dj = model->djRegion();
    double *dual = model->dualRowSolution();
    double *cost = model->costRegion();
    ClpSimplex::Status iStatus;
    const int *columnLength = matrix_->getVectorLengths();
    const CoinBigIndex *columnStart = matrix_->getVectorStarts();
    const int *row = matrix_->getIndices();
    const double *elementByColumn = matrix_->getElements();
    double infeasibilityCost = model->infeasibilityCost();
    sumDualInfeasibilities_ = 0.0;
    numberDualInfeasibilities_ = 0;
    double dualTolerance = model->dualTolerance();
    double relaxedTolerance = dualTolerance;
    // we can't really trust infeasibilities if there is dual error
    double error = std::min(1.0e-2, model->largestDualError());
    // allow tolerance at least slightly bigger than standard
    relaxedTolerance = relaxedTolerance + error;
    // but we will be using difference
    relaxedTolerance -= dualTolerance;
    sumOfRelaxedDualInfeasibilities_ = 0.0;
    for (i = 0; i < numberSets_; i++) {
      int kColumn = keyVariable_[i];
      if (kColumn < numberColumns) {
        // dj without set
        double value = cost[kColumn];
        for (CoinBigIndex j = columnStart[kColumn];
             j < columnStart[kColumn] + columnLength[kColumn]; j++) {
          int iRow = row[j];
          value -= dual[iRow] * elementByColumn[j];
        }
        // Now subtract out from all
        dj[kColumn] -= value;
        int stop = -(kColumn + 1);
        kColumn = next_[kColumn];
        while (kColumn != stop) {
          if (kColumn < 0)
            kColumn = -kColumn - 1;
          double djValue = dj[kColumn] - value;
          dj[kColumn] = djValue;
          ;
          double infeasibility = 0.0;
          iStatus = model->getStatus(kColumn);
          if (iStatus == ClpSimplex::atLowerBound) {
            if (djValue < -dualTolerance)
              infeasibility = -djValue - dualTolerance;
          } else if (iStatus == ClpSimplex::atUpperBound) {
            // at upper bound
            if (djValue > dualTolerance)
              infeasibility = djValue - dualTolerance;
          }
          if (infeasibility > 0.0) {
            sumDualInfeasibilities_ += infeasibility;
            if (infeasibility > relaxedTolerance)
              sumOfRelaxedDualInfeasibilities_ += infeasibility;
            numberDualInfeasibilities_++;
          }
          kColumn = next_[kColumn];
        }
        // check slack
        iStatus = getStatus(i);
        assert(iStatus != ClpSimplex::basic);
        double infeasibility = 0.0;
        // dj of slack is -(-1.0)value
        if (iStatus == ClpSimplex::atLowerBound) {
          if (value < -dualTolerance)
            infeasibility = -value - dualTolerance;
        } else if (iStatus == ClpSimplex::atUpperBound) {
          // at upper bound
          if (value > dualTolerance)
            infeasibility = value - dualTolerance;
        }
        if (infeasibility > 0.0) {
          sumDualInfeasibilities_ += infeasibility;
          if (infeasibility > relaxedTolerance)
            sumOfRelaxedDualInfeasibilities_ += infeasibility;
          numberDualInfeasibilities_++;
        }
      } else {
        // slack key - may not be feasible
        assert(getStatus(i) == ClpSimplex::basic);
        // negative as -1.0 for slack
        double value = -weight(i) * infeasibilityCost;
        if (value) {
          // Now subtract out from all
          int kColumn = i + numberColumns;
          int stop = -(kColumn + 1);
          kColumn = next_[kColumn];
          while (kColumn != stop) {
            if (kColumn < 0)
              kColumn = -kColumn - 1;
            double djValue = dj[kColumn] - value;
            dj[kColumn] = djValue;
            ;
            double infeasibility = 0.0;
            iStatus = model->getStatus(kColumn);
            if (iStatus == ClpSimplex::atLowerBound) {
              if (djValue < -dualTolerance)
                infeasibility = -djValue - dualTolerance;
            } else if (iStatus == ClpSimplex::atUpperBound) {
              // at upper bound
              if (djValue > dualTolerance)
                infeasibility = djValue - dualTolerance;
            }
            if (infeasibility > 0.0) {
              sumDualInfeasibilities_ += infeasibility;
              if (infeasibility > relaxedTolerance)
                sumOfRelaxedDualInfeasibilities_ += infeasibility;
              numberDualInfeasibilities_++;
            }
            kColumn = next_[kColumn];
          }
        }
      }
    }
    // and get statistics for column generation
    synchronize(model, 4);
    infeasibilityWeight_ = -1.0;
  } break;
  // Report on infeasibilities of key variables
  case 3: {
    model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_);
    model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_);
    model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_);
  } break;
  // modify costs before transposeUpdate for partial pricing
  case 4: {
    // First compute new costs etc for interesting gubs
    int iLook = 0;
    int iSet = fromIndex_[0];
    double primalTolerance = model->primalTolerance();
    const double *cost = model->costRegion();
    double *solution = model->solutionRegion();
    double infeasibilityCost = model->infeasibilityCost();
    int numberColumns = model->numberColumns();
    int numberChanged = 0;
    int *pivotVariable = model->pivotVariable();
    while (iSet >= 0) {
      int key = keyVariable_[iSet];
      double value = 0.0;
      // sum over all except key
      if ((gubType_ & 8) != 0) {
        int iColumn = next_[key];
        // sum all non-key variables
        while (iColumn >= 0) {
          value += solution[iColumn];
          iColumn = next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(key + 1);
        int iColumn = next_[key];
        // sum all non-key variables
        while (iColumn != stop) {
          if (iColumn < 0)
            iColumn = -iColumn - 1;
          value += solution[iColumn];
          iColumn = next_[iColumn];
        }
      }
      double costChange;
      double oldCost = changeCost_[iLook];
      if (key < numberColumns) {
        assert(getStatus(iSet) != ClpSimplex::basic);
        double sol;
        if (getStatus(iSet) == ClpSimplex::atUpperBound)
          sol = upper_[iSet] - value;
        else
          sol = lower_[iSet] - value;
        solution[key] = sol;
        // fix up cost
        model->nonLinearCost()->setOne(key, sol);
#ifdef CLP_DEBUG_PRINT
        printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
          cost[key], oldCost);
#endif
        costChange = cost[key] - oldCost;
      } else {
        // slack is key
        if (value > upper_[iSet] + primalTolerance) {
          setAbove(iSet);
        } else if (value < lower_[iSet] - primalTolerance) {
          setBelow(iSet);
        } else {
          setFeasible(iSet);
        }
        // negative as -1.0 for slack
        costChange = -weight(iSet) * infeasibilityCost - oldCost;
#ifdef CLP_DEBUG_PRINT
        printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
          weight(iSet) * infeasibilityCost, oldCost);
#endif
      }
      if (costChange) {
        fromIndex_[numberChanged] = iSet;
        toIndex_[iSet] = numberChanged;
        changeCost_[numberChanged++] = costChange;
      }
      iSet = fromIndex_[++iLook];
    }
    if (numberChanged || possiblePivotKey_ >= 0) {
      // first do those in list already
      int number = array->getNumElements();
      array->setPacked();
      int i;
      double *work = array->denseVector();
      int *which = array->getIndices();
      for (i = 0; i < number; i++) {
        int iRow = which[i];
        int iPivot = pivotVariable[iRow];
        if (iPivot < numberColumns) {
          int iSet = backward_[iPivot];
          if (iSet >= 0 && toIndex_[iSet] >= 0) {
            double newValue = work[i] + changeCost_[toIndex_[iSet]];
            if (!newValue)
              newValue = 1.0e-100;
            work[i] = newValue;
            // mark as done
            backward_[iPivot] = -1;
          }
        }
        if (possiblePivotKey_ == iRow) {
          double newValue = work[i] - model->dualIn();
          if (!newValue)
            newValue = 1.0e-100;
          work[i] = newValue;
          possiblePivotKey_ = -1;
        }
      }
      // now do rest and clean up
      for (i = 0; i < numberChanged; i++) {
        int iSet = fromIndex_[i];
        int key = keyVariable_[iSet];
        int iColumn = next_[key];
        double change = changeCost_[i];
        while (iColumn >= 0) {
          if (backward_[iColumn] >= 0) {
            int iRow = backToPivotRow_[iColumn];
            assert(iRow >= 0);
            work[number] = change;
            if (possiblePivotKey_ == iRow) {
              double newValue = work[number] - model->dualIn();
              if (!newValue)
                newValue = 1.0e-100;
              work[number] = newValue;
              possiblePivotKey_ = -1;
            }
            which[number++] = iRow;
          } else {
            // reset
            backward_[iColumn] = iSet;
          }
          iColumn = next_[iColumn];
        }
        toIndex_[iSet] = -1;
      }
      if (possiblePivotKey_ >= 0) {
        work[number] = -model->dualIn();
        which[number++] = possiblePivotKey_;
        possiblePivotKey_ = -1;
      }
      fromIndex_[0] = -1;
      array->setNumElements(number);
    }
  } break;
  }
}
// This is local to Gub to allow synchronization when status is good
int ClpGubMatrix::synchronize(ClpSimplex *, int)
{
  return 0;
}
/*
     general utility function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
int ClpGubMatrix::generalExpanded(ClpSimplex *model, int mode, int &number)
{
  int returnCode = 0;
  int numberColumns = model->numberColumns();
  switch (mode) {
    // Fill in pivotVariable but not for key variables
  case 0: {
    if (!next_) {
      // do ordering
      assert(!rhsOffset_);
      // create and do gub crash
      useEffectiveRhs(model, false);
    }
    int i;
    int numberBasic = number;
    // Use different array so can build from true pivotVariable_
    //int * pivotVariable = model->pivotVariable();
    int *pivotVariable = model->rowArray(0)->getIndices();
    for (i = 0; i < numberColumns; i++) {
      if (model->getColumnStatus(i) == ClpSimplex::basic) {
        int iSet = backward_[i];
        if (iSet < 0 || i != keyVariable_[iSet])
          pivotVariable[numberBasic++] = i;
      }
    }
    number = numberBasic;
    if (model->numberIterations())
      assert(number == model->numberRows());
  } break;
  // Make all key variables basic
  case 1: {
    int i;
    for (i = 0; i < numberSets_; i++) {
      int iColumn = keyVariable_[i];
      if (iColumn < numberColumns)
        model->setColumnStatus(iColumn, ClpSimplex::basic);
    }
  } break;
  // Do initial extra rows + maximum basic
  case 2: {
    returnCode = getNumRows() + 1;
    number = model->numberRows() + numberSets_;
  } break;
  // Before normal replaceColumn
  case 3: {
    int sequenceIn = model->sequenceIn();
    int sequenceOut = model->sequenceOut();
    int numberColumns = model->numberColumns();
    int numberRows = model->numberRows();
    int pivotRow = model->pivotRow();
    if (gubSlackIn_ >= 0)
      assert(sequenceIn > numberRows + numberColumns);
    if (sequenceIn == sequenceOut)
      return -1;
    int iSetIn = -1;
    int iSetOut = -1;
    if (sequenceOut < numberColumns) {
      iSetOut = backward_[sequenceOut];
    } else if (sequenceOut >= numberRows + numberColumns) {
      assert(pivotRow >= numberRows);
      int iExtra = pivotRow - numberRows;
      assert(iExtra >= 0);
      if (iSetOut < 0)
        iSetOut = fromIndex_[iExtra];
      else
        assert(iSetOut == fromIndex_[iExtra]);
    }
    if (sequenceIn < numberColumns) {
      iSetIn = backward_[sequenceIn];
    } else if (gubSlackIn_ >= 0) {
      iSetIn = gubSlackIn_;
    }
    possiblePivotKey_ = -1;
    number = 0; // say do ordinary
    int *pivotVariable = model->pivotVariable();
    if (pivotRow >= numberRows) {
      int iExtra = pivotRow - numberRows;
      //const int * length = matrix_->getVectorLengths();

      assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
      int incomingColumn = sequenceIn; // to be used in updates
      if (iSetIn != iSetOut) {
        // We need to find a possible pivot for incoming
        // look through rowArray_[1]
        int n = model->rowArray(1)->getNumElements();
        int *which = model->rowArray(1)->getIndices();
        double *array = model->rowArray(1)->denseVector();
        double bestAlpha = 1.0e-5;
        //int shortest=numberRows+1;
        for (int i = 0; i < n; i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
            if (fabs(array[i]) > fabs(bestAlpha)) {
              bestAlpha = array[i];
              possiblePivotKey_ = iRow;
            }
          }
        }
        assert(possiblePivotKey_ >= 0); // could set returnCode=4
        number = 1;
        if (sequenceIn >= numberRows + numberColumns) {
          number = 3;
          // need swap as gub slack in and must become key
          // is this best way
          int key = keyVariable_[iSetIn];
          assert(key < numberColumns);
          // check other basic
          int iColumn = next_[key];
          // set new key to be used by unpack
          keyVariable_[iSetIn] = iSetIn + numberColumns;
          // change cost in changeCost
          {
            int iLook = 0;
            int iSet = fromIndex_[0];
            while (iSet >= 0) {
              if (iSet == iSetIn) {
                changeCost_[iLook] = 0.0;
                break;
              }
              iSet = fromIndex_[++iLook];
            }
          }
          while (iColumn >= 0) {
            if (iColumn != sequenceOut) {
              // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
              printf("TTTTTry 5\n");
#endif
              int iRow = backToPivotRow_[iColumn];
              assert(iRow >= 0);
              unpack(model, model->rowArray(3), iColumn);
              model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
              double alpha = model->rowArray(3)->denseVector()[iRow];
              //if (!alpha)
              //printf("zero alpha a\n");
              int updateStatus = model->factorization()->replaceColumn(model,
                model->rowArray(2),
                model->rowArray(3),
                iRow, alpha);
              returnCode = std::max(updateStatus, returnCode);
              model->rowArray(3)->clear();
              if (returnCode)
                break;
            }
            iColumn = next_[iColumn];
          }
          if (!returnCode) {
            // now factorization looks as if key is out
            // pivot back in
#ifdef CLP_DEBUG_PRINT
            printf("TTTTTry 6\n");
#endif
            unpack(model, model->rowArray(3), key);
            model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
            pivotRow = possiblePivotKey_;
            double alpha = model->rowArray(3)->denseVector()[pivotRow];
            //if (!alpha)
            //printf("zero alpha b\n");
            int updateStatus = model->factorization()->replaceColumn(model,
              model->rowArray(2),
              model->rowArray(3),
              pivotRow, alpha);
            returnCode = std::max(updateStatus, returnCode);
            model->rowArray(3)->clear();
          }
          // restore key
          keyVariable_[iSetIn] = key;
          // now alternate column can replace key on out
          incomingColumn = pivotVariable[possiblePivotKey_];
        } else {
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTTry 4 %d\n", possiblePivotKey_);
#endif
          int updateStatus = model->factorization()->replaceColumn(model,
            model->rowArray(2),
            model->rowArray(1),
            possiblePivotKey_,
            bestAlpha);
          returnCode = std::max(updateStatus, returnCode);
          incomingColumn = pivotVariable[possiblePivotKey_];
        }

        //returnCode=4; // need swap
      } else {
        // key swap
        number = -1;
      }
      int key = keyVariable_[iSetOut];
      if (key < numberColumns)
        assert(key == sequenceOut);
      // check if any other basic
      int iColumn = next_[key];
      if (returnCode)
        iColumn = -1; // skip if error on previous
      // set new key to be used by unpack
      if (incomingColumn < numberColumns)
        keyVariable_[iSetOut] = incomingColumn;
      else
        keyVariable_[iSetOut] = iSetIn + numberColumns;
      double *cost = model->costRegion();
      if (possiblePivotKey_ < 0) {
        double dj = model->djRegion()[sequenceIn] - cost[sequenceIn];
        changeCost_[iExtra] = -dj;
#ifdef CLP_DEBUG_PRINT
        printf("modifying changeCost %d by %g - cost %g\n", iExtra, dj, cost[sequenceIn]);
#endif
      }
      while (iColumn >= 0) {
        if (iColumn != incomingColumn) {
          number = -2;
          // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTTry 1\n");
#endif
          int iRow = backToPivotRow_[iColumn];
          assert(iRow >= 0 && iRow < numberRows);
          unpack(model, model->rowArray(3), iColumn);
          model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
          double *array = model->rowArray(3)->denseVector();
          double alpha = array[iRow];
          //if (!alpha)
          //printf("zero alpha d\n");
          int updateStatus = model->factorization()->replaceColumn(model,
            model->rowArray(2),
            model->rowArray(3),
            iRow, alpha);
          returnCode = std::max(updateStatus, returnCode);
          model->rowArray(3)->clear();
          if (returnCode)
            break;
        }
        iColumn = next_[iColumn];
      }
      // restore key
      keyVariable_[iSetOut] = key;
    } else if (sequenceIn >= numberRows + numberColumns) {
      number = 2;
      //returnCode=4;
      // need swap as gub slack in and must become key
      // is this best way
      int key = keyVariable_[iSetIn];
      assert(key < numberColumns);
      // check other basic
      int iColumn = next_[key];
      // set new key to be used by unpack
      keyVariable_[iSetIn] = iSetIn + numberColumns;
      // change cost in changeCost
      {
        int iLook = 0;
        int iSet = fromIndex_[0];
        while (iSet >= 0) {
          if (iSet == iSetIn) {
            changeCost_[iLook] = 0.0;
            break;
          }
          iSet = fromIndex_[++iLook];
        }
      }
      while (iColumn >= 0) {
        if (iColumn != sequenceOut) {
          // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTry 2\n");
#endif
          int iRow = backToPivotRow_[iColumn];
          assert(iRow >= 0);
          unpack(model, model->rowArray(3), iColumn);
          model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
          double alpha = model->rowArray(3)->denseVector()[iRow];
          //if (!alpha)
          //printf("zero alpha e\n");
          int updateStatus = model->factorization()->replaceColumn(model,
            model->rowArray(2),
            model->rowArray(3),
            iRow, alpha);
          returnCode = std::max(updateStatus, returnCode);
          model->rowArray(3)->clear();
          if (returnCode)
            break;
        }
        iColumn = next_[iColumn];
      }
      if (!returnCode) {
        // now factorization looks as if key is out
        // pivot back in
#ifdef CLP_DEBUG_PRINT
        printf("TTTTTry 3\n");
#endif
        unpack(model, model->rowArray(3), key);
        model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
        double alpha = model->rowArray(3)->denseVector()[pivotRow];
        //if (!alpha)
        //printf("zero alpha f\n");
        int updateStatus = model->factorization()->replaceColumn(model,
          model->rowArray(2),
          model->rowArray(3),
          pivotRow, alpha);
        returnCode = std::max(updateStatus, returnCode);
        model->rowArray(3)->clear();
      }
      // restore key
      keyVariable_[iSetIn] = key;
    } else {
      // normal - but might as well do here
      returnCode = model->factorization()->replaceColumn(model,
        model->rowArray(2),
        model->rowArray(1),
        model->pivotRow(),
        model->alpha());
    }
  }
#ifdef CLP_DEBUG_PRINT
    printf("Update type after %d - status %d - pivot row %d\n",
      number, returnCode, model->pivotRow());
#endif
    // see if column generation says time to re-factorize
    returnCode = std::max(returnCode, synchronize(model, 5));
    number = -1; // say no need for normal replaceColumn
    break;
  // To see if can dual or primal
  case 4: {
    returnCode = 1;
  } break;
  // save status
  case 5: {
    synchronize(model, 0);
    CoinMemcpyN(status_, numberSets_, saveStatus_);
    CoinMemcpyN(keyVariable_, numberSets_, savedKeyVariable_);
  } break;
  // restore status
  case 6: {
    CoinMemcpyN(saveStatus_, numberSets_, status_);
    CoinMemcpyN(savedKeyVariable_, numberSets_, keyVariable_);
    // restore firstAvailable_
    synchronize(model, 7);
    // redo next_
    int i;
    int *last = new int[numberSets_];
    for (i = 0; i < numberSets_; i++) {
      int iKey = keyVariable_[i];
      assert(iKey >= numberColumns || backward_[iKey] == i);
      last[i] = iKey;
      // make sure basic
      //if (iKey<numberColumns)
      //model->setStatus(iKey,ClpSimplex::basic);
    }
    for (i = 0; i < numberColumns; i++) {
      int iSet = backward_[i];
      if (iSet >= 0) {
        next_[last[iSet]] = i;
        last[iSet] = i;
      }
    }
    for (i = 0; i < numberSets_; i++) {
      next_[last[i]] = -(keyVariable_[i] + 1);
      redoSet(model, keyVariable_[i], keyVariable_[i], i);
    }
    delete[] last;
    // redo pivotVariable
    int *pivotVariable = model->pivotVariable();
    int iRow;
#ifndef NDEBUG
    int numberBasic = 0;
#endif
    int numberRows = model->numberRows();
    for (iRow = 0; iRow < numberRows; iRow++) {
      if (model->getRowStatus(iRow) == ClpSimplex::basic) {
#ifndef NDEBUG
        numberBasic++;
#endif
        pivotVariable[iRow] = iRow + numberColumns;
      } else {
        pivotVariable[iRow] = -1;
      }
    }
    i = 0;
    int iColumn;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (model->getStatus(iColumn) == ClpSimplex::basic) {
        int iSet = backward_[iColumn];
        if (iSet < 0 || keyVariable_[iSet] != iColumn) {
          while (pivotVariable[i] >= 0) {
            i++;
            assert(i < numberRows);
          }
          pivotVariable[i] = iColumn;
          backToPivotRow_[iColumn] = i;
#ifndef NDEBUG
          numberBasic++;
#endif
        }
      }
    }
    assert(numberBasic == numberRows);
    rhsOffset(model, true);
  } break;
  // flag a variable
  case 7: {
    assert(number == model->sequenceIn());
    synchronize(model, 1);
    synchronize(model, 8);
  } break;
  // unflag all variables
  case 8: {
    returnCode = synchronize(model, 2);
  } break;
  // redo costs in primal
  case 9: {
    returnCode = synchronize(model, 3);
  } break;
  // return 1 if there may be changing bounds on variable (column generation)
  case 10: {
    returnCode = synchronize(model, 6);
  } break;
  // make sure set is clean
  case 11: {
    assert(number == model->sequenceIn());
    returnCode = synchronize(model, 8);
  } break;
  default:
    break;
  }
  return returnCode;
}
// Sets up an effective RHS and does gub crash if needed
void ClpGubMatrix::useEffectiveRhs(ClpSimplex *model, bool cheapest)
{
  // Do basis - cheapest or slack if feasible (unless cheapest set)
  int longestSet = 0;
  int iSet;
  for (iSet = 0; iSet < numberSets_; iSet++)
    longestSet = std::max(longestSet, end_[iSet] - start_[iSet]);

  double *upper = new double[longestSet + 1];
  double *cost = new double[longestSet + 1];
  double *lower = new double[longestSet + 1];
  double *solution = new double[longestSet + 1];
  assert(!next_);
  int numberColumns = getNumCols();
  const int *columnLength = matrix_->getVectorLengths();
  const double *columnLower = model->lowerRegion();
  const double *columnUpper = model->upperRegion();
  double *columnSolution = model->solutionRegion();
  const double *objective = model->costRegion();
  int numberRows = getNumRows();
  toIndex_ = new int[numberSets_];
  for (iSet = 0; iSet < numberSets_; iSet++)
    toIndex_[iSet] = -1;
  fromIndex_ = new int[getNumRows() + 1];
  double tolerance = model->primalTolerance();
  bool noNormalBounds = true;
  gubType_ &= ~8;
  bool gotBasis = false;
  for (iSet = 0; iSet < numberSets_; iSet++) {
    if (keyVariable_[iSet] < numberColumns)
      gotBasis = true;
    CoinBigIndex j;
    CoinBigIndex iStart = start_[iSet];
    CoinBigIndex iEnd = end_[iSet];
    for (j = iStart; j < iEnd; j++) {
      if (columnLower[j] && columnLower[j] > -1.0e20)
        noNormalBounds = false;
      if (columnUpper[j] && columnUpper[j] < 1.0e20)
        noNormalBounds = false;
    }
  }
  if (noNormalBounds)
    gubType_ |= 8;
  if (!gotBasis) {
    for (iSet = 0; iSet < numberSets_; iSet++) {
      CoinBigIndex j;
      int numberBasic = 0;
      int iBasic = -1;
      CoinBigIndex iStart = start_[iSet];
      CoinBigIndex iEnd = end_[iSet];
      // find one with smallest length
      int smallest = numberRows + 1;
      double value = 0.0;
      for (j = iStart; j < iEnd; j++) {
        if (model->getStatus(j) == ClpSimplex::basic) {
          if (columnLength[j] < smallest) {
            smallest = columnLength[j];
            iBasic = j;
          }
          numberBasic++;
        }
        value += columnSolution[j];
      }
      bool done = false;
      if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
        if (getStatus(iSet) == ClpSimplex::basic)
          iBasic = iSet + numberColumns; // slack key - use
        done = true;
      } else if (numberBasic == 1) {
        // see if can be key
        double thisSolution = columnSolution[iBasic];
        if (thisSolution > columnUpper[iBasic]) {
          value -= thisSolution - columnUpper[iBasic];
          thisSolution = columnUpper[iBasic];
          columnSolution[iBasic] = thisSolution;
        }
        if (thisSolution < columnLower[iBasic]) {
          value -= thisSolution - columnLower[iBasic];
          thisSolution = columnLower[iBasic];
          columnSolution[iBasic] = thisSolution;
        }
        // try setting slack to a bound
        assert(upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
        double cost1 = COIN_DBL_MAX;
        int whichBound = -1;
        if (upper_[iSet] < 1.0e20) {
          // try slack at ub
          double newBasic = thisSolution + upper_[iSet] - value;
          if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) {
            // can go
            whichBound = 1;
            cost1 = newBasic * objective[iBasic];
            // But if exact then may be good solution
            if (fabs(upper_[iSet] - value) < tolerance)
              cost1 = -COIN_DBL_MAX;
          }
        }
        if (lower_[iSet] > -1.0e20) {
          // try slack at lb
          double newBasic = thisSolution + lower_[iSet] - value;
          if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) {
            // can go but is it cheaper
            double cost2 = newBasic * objective[iBasic];
            // But if exact then may be good solution
            if (fabs(lower_[iSet] - value) < tolerance)
              cost2 = -COIN_DBL_MAX;
            if (cost2 < cost1)
              whichBound = 0;
          }
        }
        if (whichBound != -1) {
          // key
          done = true;
          if (whichBound) {
            // slack to upper
            columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
            setStatus(iSet, ClpSimplex::atUpperBound);
          } else {
            // slack to lower
            columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
            setStatus(iSet, ClpSimplex::atLowerBound);
          }
        }
      }
      if (!done) {
        if (!cheapest) {
          // see if slack can be key
          if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
            done = true;
            setStatus(iSet, ClpSimplex::basic);
            iBasic = iSet + numberColumns;
          }
        }
        if (!done) {
          // set non basic if there was one
          if (iBasic >= 0)
            model->setStatus(iBasic, ClpSimplex::atLowerBound);
          // find cheapest
          int numberInSet = iEnd - iStart;
          CoinMemcpyN(columnLower + iStart, numberInSet, lower);
          CoinMemcpyN(columnUpper + iStart, numberInSet, upper);
          CoinMemcpyN(columnSolution + iStart, numberInSet, solution);
          // and slack
          iBasic = numberInSet;
          solution[iBasic] = -value;
          lower[iBasic] = -upper_[iSet];
          upper[iBasic] = -lower_[iSet];
          int kphase;
          if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
            // feasible
            kphase = 1;
            cost[iBasic] = 0.0;
            CoinMemcpyN(objective + iStart, numberInSet, cost);
          } else {
            // infeasible
            kphase = 0;
            // remember bounds are flipped so opposite to natural
            if (value < lower_[iSet] - tolerance)
              cost[iBasic] = 1.0;
            else
              cost[iBasic] = -1.0;
            CoinZeroN(cost, numberInSet);
          }
          double dualTolerance = model->dualTolerance();
          for (int iphase = kphase; iphase < 2; iphase++) {
            if (iphase) {
              cost[numberInSet] = 0.0;
              CoinMemcpyN(objective + iStart, numberInSet, cost);
            }
            // now do one row lp
            bool improve = true;
            while (improve) {
              improve = false;
              double dual = cost[iBasic];
              int chosen = -1;
              double best = dualTolerance;
              int way = 0;
              for (int i = 0; i <= numberInSet; i++) {
                double dj = cost[i] - dual;
                double improvement = 0.0;
                if (iphase || i < numberInSet)
                  assert(solution[i] >= lower[i] && solution[i] <= upper[i]);
                if (dj > dualTolerance)
                  improvement = dj * (solution[i] - lower[i]);
                else if (dj < -dualTolerance)
                  improvement = dj * (solution[i] - upper[i]);
                if (improvement > best) {
                  best = improvement;
                  chosen = i;
                  if (dj < 0.0) {
                    way = 1;
                  } else {
                    way = -1;
                  }
                }
              }
              if (chosen >= 0) {
                improve = true;
                // now see how far
                if (way > 0) {
                  // incoming increasing so basic decreasing
                  // if phase 0 then go to nearest bound
                  double distance = upper[chosen] - solution[chosen];
                  double basicDistance;
                  if (!iphase) {
                    assert(iBasic == numberInSet);
                    assert(solution[iBasic] > upper[iBasic]);
                    basicDistance = solution[iBasic] - upper[iBasic];
                  } else {
                    basicDistance = solution[iBasic] - lower[iBasic];
                  }
                  // need extra coding for unbounded
                  assert(std::min(distance, basicDistance) < 1.0e20);
                  if (distance > basicDistance) {
                    // incoming becomes basic
                    solution[chosen] += basicDistance;
                    if (!iphase)
                      solution[iBasic] = upper[iBasic];
                    else
                      solution[iBasic] = lower[iBasic];
                    iBasic = chosen;
                  } else {
                    // flip
                    solution[chosen] = upper[chosen];
                    solution[iBasic] -= distance;
                  }
                } else {
                  // incoming decreasing so basic increasing
                  // if phase 0 then go to nearest bound
                  double distance = solution[chosen] - lower[chosen];
                  double basicDistance;
                  if (!iphase) {
                    assert(iBasic == numberInSet);
                    assert(solution[iBasic] < lower[iBasic]);
                    basicDistance = lower[iBasic] - solution[iBasic];
                  } else {
                    basicDistance = upper[iBasic] - solution[iBasic];
                  }
                  // need extra coding for unbounded - for now just exit
                  if (std::min(distance, basicDistance) > 1.0e20) {
                    printf("unbounded on set %d\n", iSet);
                    iphase = 1;
                    iBasic = numberInSet;
                    break;
                  }
                  if (distance > basicDistance) {
                    // incoming becomes basic
                    solution[chosen] -= basicDistance;
                    if (!iphase)
                      solution[iBasic] = lower[iBasic];
                    else
                      solution[iBasic] = upper[iBasic];
                    iBasic = chosen;
                  } else {
                    // flip
                    solution[chosen] = lower[chosen];
                    solution[iBasic] += distance;
                  }
                }
                if (!iphase) {
                  if (iBasic < numberInSet)
                    break; // feasible
                  else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic])
                    break; // feasible (on flip)
                }
              }
            }
          }
          // convert iBasic back and do bounds
          if (iBasic == numberInSet) {
            // slack basic
            setStatus(iSet, ClpSimplex::basic);
            iBasic = iSet + numberColumns;
          } else {
            iBasic += start_[iSet];
            model->setStatus(iBasic, ClpSimplex::basic);
            // remember bounds flipped
            if (upper[numberInSet] == lower[numberInSet])
              setStatus(iSet, ClpSimplex::isFixed);
            else if (solution[numberInSet] == upper[numberInSet])
              setStatus(iSet, ClpSimplex::atLowerBound);
            else if (solution[numberInSet] == lower[numberInSet])
              setStatus(iSet, ClpSimplex::atUpperBound);
            else
              abort();
          }
          for (j = iStart; j < iEnd; j++) {
            if (model->getStatus(j) != ClpSimplex::basic) {
              int inSet = j - iStart;
              columnSolution[j] = solution[inSet];
              if (upper[inSet] == lower[inSet])
                model->setStatus(j, ClpSimplex::isFixed);
              else if (solution[inSet] == upper[inSet])
                model->setStatus(j, ClpSimplex::atUpperBound);
              else if (solution[inSet] == lower[inSet])
                model->setStatus(j, ClpSimplex::atLowerBound);
            }
          }
        }
      }
      keyVariable_[iSet] = iBasic;
    }
  }
  delete[] lower;
  delete[] solution;
  delete[] upper;
  delete[] cost;
  // make sure matrix is in good shape
  matrix_->orderMatrix();
  // create effective rhs
  delete[] rhsOffset_;
  rhsOffset_ = new double[numberRows];
  delete[] next_;
  next_ = new int[numberColumns + numberSets_ + 2 * longestSet];
  char *mark = new char[numberColumns];
  memset(mark, 0, numberColumns);
  for (int iColumn = 0; iColumn < numberColumns; iColumn++)
    next_[iColumn] = COIN_INT_MAX;
  int i;
  int *keys = new int[numberSets_];
  for (i = 0; i < numberSets_; i++)
    keys[i] = COIN_INT_MAX;
  // set up chains
  for (i = 0; i < numberColumns; i++) {
    if (model->getStatus(i) == ClpSimplex::basic)
      mark[i] = 1;
    int iSet = backward_[i];
    if (iSet >= 0) {
      int iNext = keys[iSet];
      next_[i] = iNext;
      keys[iSet] = i;
    }
  }
  for (i = 0; i < numberSets_; i++) {
    int j;
    if (getStatus(i) != ClpSimplex::basic) {
      // make sure fixed if it is
      if (upper_[i] == lower_[i])
        setStatus(i, ClpSimplex::isFixed);
      // slack not key - choose one with smallest length
      int smallest = numberRows + 1;
      int key = -1;
      j = keys[i];
      if (j != COIN_INT_MAX) {
        while (1) {
          if (mark[j] && columnLength[j] < smallest && !gotBasis) {
            key = j;
            smallest = columnLength[j];
          }
          if (next_[j] != COIN_INT_MAX) {
            j = next_[j];
          } else {
            // correct end
            next_[j] = -(keys[i] + 1);
            break;
          }
        }
      } else {
        next_[i + numberColumns] = -(numberColumns + i + 1);
      }
      if (gotBasis)
        key = keyVariable_[i];
      if (key >= 0) {
        keyVariable_[i] = key;
      } else {
        // nothing basic - make slack key
        //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
        // fudge to avoid const problem
        status_[i] = 1;
      }
    } else {
      // slack key
      keyVariable_[i] = numberColumns + i;
      int j;
      double sol = 0.0;
      j = keys[i];
      if (j != COIN_INT_MAX) {
        while (1) {
          sol += columnSolution[j];
          if (next_[j] != COIN_INT_MAX) {
            j = next_[j];
          } else {
            // correct end
            next_[j] = -(keys[i] + 1);
            break;
          }
        }
      } else {
        next_[i + numberColumns] = -(numberColumns + i + 1);
      }
      if (sol > upper_[i] + tolerance) {
        setAbove(i);
      } else if (sol < lower_[i] - tolerance) {
        setBelow(i);
      } else {
        setFeasible(i);
      }
    }
    // Create next_
    int key = keyVariable_[i];
    redoSet(model, key, keys[i], i);
  }
  delete[] keys;
  delete[] mark;
  rhsOffset(model, true);
}
// redoes next_ for a set.
void ClpGubMatrix::redoSet(ClpSimplex *model, int newKey, int oldKey, int iSet)
{
  int numberColumns = model->numberColumns();
  int *save = next_ + numberColumns + numberSets_;
  int number = 0;
  int stop = -(oldKey + 1);
  int j = next_[oldKey];
  while (j != stop) {
    if (j < 0)
      j = -j - 1;
    if (j != newKey)
      save[number++] = j;
    j = next_[j];
  }
  // and add oldkey
  if (newKey != oldKey)
    save[number++] = oldKey;
  // now do basic
  int lastMarker = -(newKey + 1);
  keyVariable_[iSet] = newKey;
  next_[newKey] = lastMarker;
  int last = newKey;
  for (j = 0; j < number; j++) {
    int iColumn = save[j];
    if (iColumn < numberColumns) {
      if (model->getStatus(iColumn) == ClpSimplex::basic) {
        next_[last] = iColumn;
        next_[iColumn] = lastMarker;
        last = iColumn;
      }
    }
  }
  // now add in non-basic
  for (j = 0; j < number; j++) {
    int iColumn = save[j];
    if (iColumn < numberColumns) {
      if (model->getStatus(iColumn) != ClpSimplex::basic) {
        next_[last] = -(iColumn + 1);
        next_[iColumn] = lastMarker;
        last = iColumn;
      }
    }
  }
}
/* Returns effective RHS if it is being used.  This is used for long problems
   or big gub or anywhere where going through full columns is
   expensive.  This may re-compute */
double *
ClpGubMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh, bool
#ifdef CLP_DEBUG
                                                                check
#endif
)
{
  //forceRefresh=true;
  if (rhsOffset_) {
#ifdef CLP_DEBUG
    if (check) {
      // no need - but check anyway
      // zero out basic
      int numberRows = model->numberRows();
      int numberColumns = model->numberColumns();
      double *solution = new double[numberColumns];
      double *rhs = new double[numberRows];
      CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
      CoinZeroN(rhs, numberRows);
      int iRow;
      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
          solution[iColumn] = 0.0;
      }
      for (int iSet = 0; iSet < numberSets_; iSet++) {
        int iColumn = keyVariable_[iSet];
        if (iColumn < numberColumns)
          solution[iColumn] = 0.0;
      }
      times(-1.0, solution, rhs);
      delete[] solution;
      const double *columnSolution = model->solutionRegion();
      // and now subtract out non basic
      ClpSimplex::Status iStatus;
      for (int iSet = 0; iSet < numberSets_; iSet++) {
        int iColumn = keyVariable_[iSet];
        if (iColumn < numberColumns) {
          double b = 0.0;
          // key is structural - where is slack
          iStatus = getStatus(iSet);
          assert(iStatus != ClpSimplex::basic);
          if (iStatus == ClpSimplex::atLowerBound)
            b = lower_[iSet];
          else
            b = upper_[iSet];
          // subtract out others at bounds
          if ((gubType_ & 8) == 0) {
            int stop = -(iColumn + 1);
            int jColumn = next_[iColumn];
            // sum all non-basic variables - first skip basic
            while (jColumn >= 0)
              jColumn = next_[jColumn];
            while (jColumn != stop) {
              assert(jColumn < 0);
              jColumn = -jColumn - 1;
              b -= columnSolution[jColumn];
              jColumn = next_[jColumn];
            }
          }
          // subtract out
          ClpPackedMatrix::add(model, rhs, iColumn, -b);
        }
      }
      for (iRow = 0; iRow < numberRows; iRow++) {
        if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
          printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
      }
      delete[] rhs;
    }
#endif
    if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) {
      // zero out basic
      int numberRows = model->numberRows();
      int numberColumns = model->numberColumns();
      double *solution = new double[numberColumns];
      CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
      CoinZeroN(rhsOffset_, numberRows);
      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
          solution[iColumn] = 0.0;
      }
      int iSet;
      for (iSet = 0; iSet < numberSets_; iSet++) {
        int iColumn = keyVariable_[iSet];
        if (iColumn < numberColumns)
          solution[iColumn] = 0.0;
      }
      times(-1.0, solution, rhsOffset_);
      delete[] solution;
      lastRefresh_ = model->numberIterations();
      const double *columnSolution = model->solutionRegion();
      // and now subtract out non basic
      ClpSimplex::Status iStatus;
      for (iSet = 0; iSet < numberSets_; iSet++) {
        int iColumn = keyVariable_[iSet];
        if (iColumn < numberColumns) {
          double b = 0.0;
          // key is structural - where is slack
          iStatus = getStatus(iSet);
          assert(iStatus != ClpSimplex::basic);
          if (iStatus == ClpSimplex::atLowerBound)
            b = lower_[iSet];
          else
            b = upper_[iSet];
          // subtract out others at bounds
          if ((gubType_ & 8) == 0) {
            int stop = -(iColumn + 1);
            int jColumn = next_[iColumn];
            // sum all non-basic variables - first skip basic
            while (jColumn >= 0)
              jColumn = next_[jColumn];
            while (jColumn != stop) {
              assert(jColumn < 0);
              jColumn = -jColumn - 1;
              b -= columnSolution[jColumn];
              jColumn = next_[jColumn];
            }
          }
          // subtract out
          if (b)
            ClpPackedMatrix::add(model, rhsOffset_, iColumn, -b);
        }
      }
    }
  }
  return rhsOffset_;
}
/*
   update information for a pivot (and effective rhs)
*/
int ClpGubMatrix::updatePivot(ClpSimplex *model, double oldInValue, double /*oldOutValue*/)
{
  int sequenceIn = model->sequenceIn();
  int sequenceOut = model->sequenceOut();
  double *solution = model->solutionRegion();
  int numberColumns = model->numberColumns();
  int numberRows = model->numberRows();
  int pivotRow = model->pivotRow();
  int iSetIn;
  // Correct sequence in
  trueSequenceIn_ = sequenceIn;
  if (sequenceIn < numberColumns) {
    iSetIn = backward_[sequenceIn];
  } else if (sequenceIn < numberColumns + numberRows) {
    iSetIn = -1;
  } else {
    iSetIn = gubSlackIn_;
    trueSequenceIn_ = numberColumns + numberRows + iSetIn;
  }
  int iSetOut = -1;
  trueSequenceOut_ = sequenceOut;
  if (sequenceOut < numberColumns) {
    iSetOut = backward_[sequenceOut];
  } else if (sequenceOut >= numberRows + numberColumns) {
    assert(pivotRow >= numberRows);
    int iExtra = pivotRow - numberRows;
    assert(iExtra >= 0);
    if (iSetOut < 0)
      iSetOut = fromIndex_[iExtra];
    else
      assert(iSetOut == fromIndex_[iExtra]);
    trueSequenceOut_ = numberColumns + numberRows + iSetOut;
  }
  if (rhsOffset_) {
    // update effective rhs
    if (sequenceIn == sequenceOut) {
      assert(sequenceIn < numberRows + numberColumns); // should be easy to deal with
      if (sequenceIn < numberColumns)
        add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]);
    } else {
      if (sequenceIn < numberColumns) {
        // we need to test if WILL be key
        ClpPackedMatrix::add(model, rhsOffset_, sequenceIn, oldInValue);
        if (iSetIn >= 0) {
          // old contribution to rhsOffset_
          int key = keyVariable_[iSetIn];
          if (key < numberColumns) {
            double oldB = 0.0;
            ClpSimplex::Status iStatus = getStatus(iSetIn);
            if (iStatus == ClpSimplex::atLowerBound)
              oldB = lower_[iSetIn];
            else
              oldB = upper_[iSetIn];
            // subtract out others at bounds
            if ((gubType_ & 8) == 0) {
              int stop = -(key + 1);
              int iColumn = next_[key];
              // skip basic
              while (iColumn >= 0)
                iColumn = next_[iColumn];
              // sum all non-key variables
              while (iColumn != stop) {
                assert(iColumn < 0);
                iColumn = -iColumn - 1;
                if (iColumn == sequenceIn)
                  oldB -= oldInValue;
                else if (iColumn != sequenceOut)
                  oldB -= solution[iColumn];
                iColumn = next_[iColumn];
              }
            }
            if (oldB)
              ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
          }
        }
      } else if (sequenceIn < numberRows + numberColumns) {
        //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
      } else {
#ifdef CLP_DEBUG_PRINT
        printf("** in is key slack %d\n", sequenceIn);
#endif
        // old contribution to rhsOffset_
        int key = keyVariable_[iSetIn];
        if (key < numberColumns) {
          double oldB = 0.0;
          ClpSimplex::Status iStatus = getStatus(iSetIn);
          if (iStatus == ClpSimplex::atLowerBound)
            oldB = lower_[iSetIn];
          else
            oldB = upper_[iSetIn];
          // subtract out others at bounds
          if ((gubType_ & 8) == 0) {
            int stop = -(key + 1);
            int iColumn = next_[key];
            // skip basic
            while (iColumn >= 0)
              iColumn = next_[iColumn];
            // sum all non-key variables
            while (iColumn != stop) {
              assert(iColumn < 0);
              iColumn = -iColumn - 1;
              if (iColumn != sequenceOut)
                oldB -= solution[iColumn];
              iColumn = next_[iColumn];
            }
          }
          if (oldB)
            ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
        }
      }
      if (sequenceOut < numberColumns) {
        ClpPackedMatrix::add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]);
        if (iSetOut >= 0) {
          // old contribution to rhsOffset_
          int key = keyVariable_[iSetOut];
          if (key < numberColumns && iSetIn != iSetOut) {
            double oldB = 0.0;
            ClpSimplex::Status iStatus = getStatus(iSetOut);
            if (iStatus == ClpSimplex::atLowerBound)
              oldB = lower_[iSetOut];
            else
              oldB = upper_[iSetOut];
            // subtract out others at bounds
            if ((gubType_ & 8) == 0) {
              int stop = -(key + 1);
              int iColumn = next_[key];
              // skip basic
              while (iColumn >= 0)
                iColumn = next_[iColumn];
              // sum all non-key variables
              while (iColumn != stop) {
                assert(iColumn < 0);
                iColumn = -iColumn - 1;
                if (iColumn == sequenceIn)
                  oldB -= oldInValue;
                else if (iColumn != sequenceOut)
                  oldB -= solution[iColumn];
                iColumn = next_[iColumn];
              }
            }
            if (oldB)
              ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
          }
        }
      } else if (sequenceOut < numberRows + numberColumns) {
        //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
      } else {
#ifdef CLP_DEBUG_PRINT
        printf("** out is key slack %d\n", sequenceOut);
#endif
        assert(pivotRow >= numberRows);
      }
    }
  }
  int *pivotVariable = model->pivotVariable();
  // may need to deal with key
  // Also need coding to mark/allow key slack entering
  if (pivotRow >= numberRows) {
    assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
#ifdef CLP_DEBUG_PRINT
    if (sequenceIn >= numberRows + numberColumns)
      printf("key slack %d in, set out %d\n", gubSlackIn_, iSetOut);
    printf("** danger - key out for set %d in %d (set %d)\n", iSetOut, sequenceIn,
      iSetIn);
#endif
    // if slack out mark correctly
    if (sequenceOut >= numberRows + numberColumns) {
      double value = model->valueOut();
      if (value == upper_[iSetOut]) {
        setStatus(iSetOut, ClpSimplex::atUpperBound);
      } else if (value == lower_[iSetOut]) {
        setStatus(iSetOut, ClpSimplex::atLowerBound);
      } else {
        if (fabs(value - upper_[iSetOut]) < fabs(value - lower_[iSetOut])) {
          setStatus(iSetOut, ClpSimplex::atUpperBound);
        } else {
          setStatus(iSetOut, ClpSimplex::atLowerBound);
        }
      }
      if (upper_[iSetOut] == lower_[iSetOut])
        setStatus(iSetOut, ClpSimplex::isFixed);
      setFeasible(iSetOut);
    }
    if (iSetOut == iSetIn) {
      // key swap
      int key;
      if (sequenceIn >= numberRows + numberColumns) {
        key = numberColumns + iSetIn;
        setStatus(iSetIn, ClpSimplex::basic);
      } else {
        key = sequenceIn;
      }
      redoSet(model, key, keyVariable_[iSetIn], iSetIn);
    } else {
      // key was chosen
      assert(possiblePivotKey_ >= 0 && possiblePivotKey_ < numberRows);
      int key = pivotVariable[possiblePivotKey_];
      // and set incoming here
      if (sequenceIn >= numberRows + numberColumns) {
        // slack in - so use old key
        sequenceIn = keyVariable_[iSetIn];
        model->setStatus(sequenceIn, ClpSimplex::basic);
        setStatus(iSetIn, ClpSimplex::basic);
        redoSet(model, iSetIn + numberColumns, keyVariable_[iSetIn], iSetIn);
      }
      //? do not do if iSetIn<0 ? as will be done later
      pivotVariable[possiblePivotKey_] = sequenceIn;
      if (sequenceIn < numberColumns)
        backToPivotRow_[sequenceIn] = possiblePivotKey_;
      redoSet(model, key, keyVariable_[iSetOut], iSetOut);
    }
  } else {
    if (sequenceOut < numberColumns) {
      if (iSetIn >= 0 && iSetOut == iSetIn) {
        // key not out - only problem is if slack in
        int key;
        if (sequenceIn >= numberRows + numberColumns) {
          key = numberColumns + iSetIn;
          setStatus(iSetIn, ClpSimplex::basic);
          assert(pivotRow < numberRows);
          // must swap with current key
          int key = keyVariable_[iSetIn];
          model->setStatus(key, ClpSimplex::basic);
          pivotVariable[pivotRow] = key;
          backToPivotRow_[key] = pivotRow;
        } else {
          key = keyVariable_[iSetIn];
        }
        redoSet(model, key, keyVariable_[iSetIn], iSetIn);
      } else if (iSetOut >= 0) {
        // just redo set
        int key = keyVariable_[iSetOut];
        ;
        redoSet(model, key, keyVariable_[iSetOut], iSetOut);
      }
    }
  }
  if (iSetIn >= 0 && iSetIn != iSetOut) {
    int key = keyVariable_[iSetIn];
    if (sequenceIn == numberColumns + 2 * numberRows) {
      // key slack in
      assert(pivotRow < numberRows);
      // must swap with current key
      model->setStatus(key, ClpSimplex::basic);
      pivotVariable[pivotRow] = key;
      backToPivotRow_[key] = pivotRow;
      setStatus(iSetIn, ClpSimplex::basic);
      key = iSetIn + numberColumns;
    }
    // redo set to allow for new one
    redoSet(model, key, keyVariable_[iSetIn], iSetIn);
  }
  // update pivot
  if (sequenceIn < numberColumns) {
    if (pivotRow < numberRows) {
      backToPivotRow_[sequenceIn] = pivotRow;
    } else {
      if (possiblePivotKey_ >= 0) {
        assert(possiblePivotKey_ < numberRows);
        backToPivotRow_[sequenceIn] = possiblePivotKey_;
        pivotVariable[possiblePivotKey_] = sequenceIn;
      }
    }
  } else if (sequenceIn >= numberRows + numberColumns) {
    // key in - something should have been done before
    int key = keyVariable_[iSetIn];
    assert(key == numberColumns + iSetIn);
    //pivotVariable[pivotRow]=key;
    //backToPivotRow_[key]=pivotRow;
    //model->setStatus(key,ClpSimplex::basic);
    //key=numberColumns+iSetIn;
    setStatus(iSetIn, ClpSimplex::basic);
    redoSet(model, key, keyVariable_[iSetIn], iSetIn);
  }
#ifdef CLP_DEBUG
  {
    char *xx = new char[numberColumns + numberRows];
    memset(xx, 0, numberRows + numberColumns);
    for (int i = 0; i < numberRows; i++) {
      int iPivot = pivotVariable[i];
      assert(iPivot < numberRows + numberColumns);
      assert(!xx[iPivot]);
      xx[iPivot] = 1;
      if (iPivot < numberColumns) {
        int iBack = backToPivotRow_[iPivot];
        assert(i == iBack);
      }
    }
    delete[] xx;
  }
#endif
  if (rhsOffset_) {
    // update effective rhs
    if (sequenceIn != sequenceOut) {
      if (sequenceIn < numberColumns) {
        if (iSetIn >= 0) {
          // new contribution to rhsOffset_
          int key = keyVariable_[iSetIn];
          if (key < numberColumns) {
            double newB = 0.0;
            ClpSimplex::Status iStatus = getStatus(iSetIn);
            if (iStatus == ClpSimplex::atLowerBound)
              newB = lower_[iSetIn];
            else
              newB = upper_[iSetIn];
            // subtract out others at bounds
            if ((gubType_ & 8) == 0) {
              int stop = -(key + 1);
              int iColumn = next_[key];
              // skip basic
              while (iColumn >= 0)
                iColumn = next_[iColumn];
              // sum all non-key variables
              while (iColumn != stop) {
                assert(iColumn < 0);
                iColumn = -iColumn - 1;
                newB -= solution[iColumn];
                iColumn = next_[iColumn];
              }
            }
            if (newB)
              ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
          }
        }
      }
      if (iSetOut >= 0) {
        // new contribution to rhsOffset_
        int key = keyVariable_[iSetOut];
        if (key < numberColumns && iSetIn != iSetOut) {
          double newB = 0.0;
          ClpSimplex::Status iStatus = getStatus(iSetOut);
          if (iStatus == ClpSimplex::atLowerBound)
            newB = lower_[iSetOut];
          else
            newB = upper_[iSetOut];
          // subtract out others at bounds
          if ((gubType_ & 8) == 0) {
            int stop = -(key + 1);
            int iColumn = next_[key];
            // skip basic
            while (iColumn >= 0)
              iColumn = next_[iColumn];
            // sum all non-key variables
            while (iColumn != stop) {
              assert(iColumn < 0);
              iColumn = -iColumn - 1;
              newB -= solution[iColumn];
              iColumn = next_[iColumn];
            }
          }
          if (newB)
            ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
        }
      }
    }
  }
#ifdef CLP_DEBUG
  // debug
  {
    int i;
    char *xxxx = new char[numberColumns];
    memset(xxxx, 0, numberColumns);
    for (i = 0; i < numberRows; i++) {
      int iPivot = pivotVariable[i];
      assert(model->getStatus(iPivot) == ClpSimplex::basic);
      if (iPivot < numberColumns && backward_[iPivot] >= 0)
        xxxx[iPivot] = 1;
    }
    double primalTolerance = model->primalTolerance();
    for (i = 0; i < numberSets_; i++) {
      int key = keyVariable_[i];
      double value = 0.0;
      // sum over all except key
      int iColumn = next_[key];
      // sum all non-key variables
      int k = 0;
      int stop = -(key + 1);
      while (iColumn != stop) {
        if (iColumn < 0)
          iColumn = -iColumn - 1;
        value += solution[iColumn];
        k++;
        assert(k < 100);
        assert(backward_[iColumn] == i);
        iColumn = next_[iColumn];
      }
      iColumn = next_[key];
      if (key < numberColumns) {
        // feasibility will be done later
        assert(getStatus(i) != ClpSimplex::basic);
        double sol;
        if (getStatus(i) == ClpSimplex::atUpperBound)
          sol = upper_[i] - value;
        else
          sol = lower_[i] - value;
        //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
        //     cost[key]);
        //if (fabs(sol-solution[key])>1.0e-3)
        //printf("** stored value was %g\n",solution[key]);
      } else {
        // slack is key
        double infeasibility = 0.0;
        if (value > upper_[i] + primalTolerance) {
          infeasibility = value - upper_[i] - primalTolerance;
          //setAbove(i);
        } else if (value < lower_[i] - primalTolerance) {
          infeasibility = lower_[i] - value - primalTolerance;
          //setBelow(i);
        } else {
          //setFeasible(i);
        }
        //printf("xx Value of key slack for set %d is %g\n",i,value);
      }
      while (iColumn >= 0) {
        assert(xxxx[iColumn]);
        xxxx[iColumn] = 0;
        iColumn = next_[iColumn];
      }
    }
    for (i = 0; i < numberColumns; i++) {
      if (i < numberColumns && backward_[i] >= 0) {
        assert(!xxxx[i] || i == keyVariable_[backward_[i]]);
      }
    }
    delete[] xxxx;
  }
#endif
  return 0;
}
// Switches off dj checking each factorization (for BIG models)
void ClpGubMatrix::switchOffCheck()
{
  noCheck_ = 0;
  infeasibilityWeight_ = 0.0;
}
// Correct sequence in and out to give true value
void ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/, int &sequenceIn, int &sequenceOut)
{
  if (sequenceIn != -999) {
    sequenceIn = trueSequenceIn_;
    sequenceOut = trueSequenceOut_;
  }
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/
